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ABSTRACT 

We explore the evolution of radiatively inefficient accretion disks in which nuclear reactions are dynamically 
important ('Nuclear Dominated Accretion Flows', or NuDAFs). Examples of such disks are those generated 
by the merger of a white dwarf with a neutron star or black hole, or by the collapse of a rotating star. Here 
we present two-dimensional hydrodynamic simulations that systematically explore the effect of adding a single 
nuclear reaction to a viscous torus. The equation of state, anomalous shear stress, and nuclear reactions are 
given parametric forms. Our results point to the existence of two qualitatively different regimes of NuDAF 
evolution: (1) steady accretion with quiescent burning; or (2) detonation of the disk. These outcomes are 
controlled primarily by the ratio $ of the nuclear energy released to the enthalpy at the burning radius. Disks 
detonate if ^ exceeds a critical value $ cr i t ~ 1, and if burning occurs in regions where neutrino cooling is 
unimportant. Thermonuclear runaways are seeded by the turbulent mixing of hot ash with cold fuel at the 
burning front. Disks with ^ < ^ crit do not explode, but instead power a persistent collimated outflow of 
unbound material composed primarily of ash, with a mass-loss rate that increases with 'J. We discuss the 
implications of our results for supernova-like counterparts from astrophysical events in the NuDAF regime. In 
particular, detonations following a white dwarf - neutron star merger could account for some subluminous Type 
la supernovae, such as the class defined by SN 2002cx. 

Subject headings: accretion, accretion disks — nuclear reactions, nucleosynthesis, abundances — white dwarfs 
— hydrodynamics — stars: winds, outflows — supernovae: general 



1. INTRODUCTION 

Many violent and visually spectacular astrophysical events 
involve stellar material collapsing at high rates onto a cen- 
tral compact object. Examples include the core collapse of a 
massive star to form a neutron star (NS) or black hole (BH) 
remnant (Woosley et al. 2002); the accretion-induced collapse 
of a white dwarf (WD) (Nomoto & Kondo 1991); the merger 
of binary compact objects, including various combinations of 
WDs, NSs, and BHs; the binary merger of a Helium star with 
a NS or BH (e.g., Fryer & Woosley 1998); and the inspiral 
and tidal disruption of a planet that is engulfed by the enve- 
lope of its host star (e.g., Siess & Livio 1999). Although only 
some of these events have yet been unambiguously observed, 
this situation may change soon as the result of new transient 
surveys coming online across the electromagnetic spectrum. 

In many of the above examples, angular momentum places 
a key dynamical role, such that a significant fraction of the 
infalling matter forms a rotationally-supported torus around 
the central object. If the accretion rate through such a disk 
is sufficiently high, then the heat generated is trapped in the 
flow instead of being efficiently radiated away (e.g., Chevalier 
1993). Such disks are examples of what is generally termed a 
Radiatively Inefficient Accretion Flow (RIAF). 4 

Previous theoretical work on RIAFs has focused on under- 
standing their spatial structure and on determining the pro- 
cesses that control the radial transport of mass, energy, and 
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angular momentum through the disk. The ADAF model 
(Narayan & Yi 1994, 1995) assumes that angular momen- 
tum is transported outward by turbulent stresses (mediated by, 
e.g., the magnetorotational instability [MRI]; Balbus & Haw- 
ley 1998) and that thermal energy is advected through the in- 
ner boundary of the disk at small radii. The CDAF model 
(Quataert & Gruzinov 2000) postulates that outward angular 
momentum transport is offset by inward transport by convec- 
tion. The resulting accretion rate is much lower than that aris- 
ing from ADAFs, though all of the disk material is eventually 
accreted by the central object. In contrast, the ADIOS model 
(Blandford & Begelman 1999) assumes that most of the ac- 
cretion energy powers an unbound outflow from the disk, such 
that the majority of the material never reaches the central ob- 
ject. 

Numerous two- and three dimensional global simulations 
of RIAFs have been performed in hydrodynamics (2D - Stone 
et al. 1999, Igumenshchev & Abramowicz 2000; 3D - Igu- 
menshchev et al. 2000) as well as in magnetohydrodynamics 
[MHD] (2D - Stone & Pringle 2001; 3D - Hawley et al. 2001, 
Machida et al. 2001, Hawley & Balbus 2002, Igumenshchev 
et al. 2003, Pen et al. 2003, Pang et al. 201 1, McKinney et al. 
2012, Narayan et al. 2012). A key conclusion of most such 
simulations to date is that the net accretion rate through the 
inner edge of the torus at r = r m is substantially reduced (by a 
factor ~ r m /Ro) from the feeding rate set by the outer torus at 
r = Rq, consistent with both CDAF and ADIOS models (e.g., 
Yuan et al. 2012). However, the relative importance of truly 
unbound outflows versus simply large-scale bound convective 
motions remains under debate (e.g., Abramowicz et al. 2000). 

Here we focus on another key feature of RIAFs in collapsar 
and merging compact object scenarios: the nuclear composi- 
tion of the accreting material and its effect on the dynamics 
of the accretion disk. If the potential well of the primary has 
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an appropriate depth, then the inflowing matter becomes suffi- 
ciently hot and dense for nuclear reactions to become dynami- 
cally important and to generate radioactive material. Depend- 
ing on how nuclear burning affects the accretion rate onto the 
central object, or the properties of outflows from the disk, this 
can greatly impact predictions for the observational signatures 
of these events. 

Metzger (2012; hereafter M12) constructed a steady-state, 
height-integrated model of accretion following the tidal dis- 
ruption of a WD by a BH or NS, including the effects of 
nuclear burning on both the thermodynamics and composi- 
tion of the accreting material. He adopted an ADIOS-type 
model, in which the properties of the outflow from each ra- 
dius in the disk are calculated by locally balancing wind cool- 
ing with viscous and nuclear heating. In this model, the addi- 
tional heating from nuclear burning acts to enhance the mass 
loss rate in the outflow over its nominal value without nuclear 
burning. Ml 2 showed that at radii r < 10 9 cm, the midplane 
density and temperature become sufficiently high to burn the 
initial white dwarf material into increasingly heavier elements 
(e.g., Mg, Si, S, Ca, Fe, and Ni) at sequentially smaller radii. 
The disk structure thus resembles the onion-like structure of 
a massive star, but one which evolves on a much more rapid 
timescale. The combined outflow from all annuli in the disk 
is composed primarily of unburnt WD material from the outer 
disk, along with a smaller fraction of intermediate mass ele- 
ments and radioactive 56 Ni originating from smaller radii. 

In some regions of the disk Ml 2 found that the energy re- 
leased by nuclear reactions is at least comparable to that re- 
leased viscously, thus motivating him to term this novel accre- 
tion regime a Nuclear Dominated Accretion Flow (NuDAF). 
Because RIAFs are already only marginally bound, even a 
modest additional source of energy can have a significant 
impact on its structure. In other words, the properties of 
NuDAFs could in principle differ markedly from the predic- 
tions of normal ADAF/CDAF/ADIOS models. 

Given its many simplifying assumptions, the M12 model 
cannot address several key questions relevant to NuDAFs. 
These include (1) the true effects of nuclear burning on 
the quantity and composition of disk outflows; (2) multi- 
dimensional effects, such as radial convective mixing or the 
diffusion of burnt fuel upstream; (3) the long-term, global 
evolution of the disk; and (4) the thermal stability of the 
inflow, which is subtle due to the sensitive temperature de- 
pendence of the nuclear reaction rates. M12 found that ther- 
mal stability depends critically on what is assumed about the 
pressure dependence of wind cooling. If the disk becomes 
thermally unstable, possible outcomes range from a complex 
'limit cycle' behavior (e.g., periods of inflow followed by in- 
tense rapid burning) to a complete dynamical explosion. 

In this series of papers we explore the effect of nuclear 
burning on the structure and evolution of RIAFs by means 
of axisymmetric hydrodynamic simulations. In paper I 
(this work) we explore the general dynamical properties of 
NuDAFs and their outflows, and compare them to the known 
case where nuclear burning is absent (e.g., Stone et al. 1999). 
As in Ml 2, our study focuses on disks created by the tidal 
disruption of a WD by a NS or BH binary companion, in part 
because the global properties of the torus are relatively well- 
constrained by the parameters of the binary. However, many 
of the conclusions we reach should apply to other similar as- 
trophysical events, such as the collapse of a rotating star. 

In this initial study we make several simplifications in order 



to clearly identify the processes controlling the dynamics, and 
to allow an efficient exploration of parameter space. The most 
important of these approximations is to parameterize the an- 
gular momentum transport in the disk by an anomalous shear 
viscosity. We also adopt an ideal gas equation of state with a 
single adiabatic index, and for ease include only one nuclear 
reaction. In a subsequent paper we will use a realistic equa- 
tion of state (EOS) and extend our nuclear reaction chain to a 
full a-reaction network, thus enabling a more reliable predic- 
tion of observational signatures. 

The paper is organized as follows. In §2 we describe the 
physics included in our model. Details on the numerical im- 
plementation are provided in §3. Results are presented in 
§4, beginning with an overview of a fiducial model, followed 
by sub-sections on exploding and quiescent models. Impli- 
cations from our results are discussed in §5, and a bulleted 
summary of our conclusions is given in §6. Conversion of 
physical to code units and reaction rates employed are de- 
scribed in Appendix A. Appendix B documents tests of our 
numerical code. An analytic model for the conditions leading 
to disk detonation via the Zel'dovich mechanism is provided 
in Appendix C. 

Readers uninterested in technical details are encouraged to 
skip directly to §4. 1 . 

2. PHYSICAL MODEL 

As a clear physical testbed to study NuDAFs, we focus on 
disks that form when a WD is tidally disrupted by a compan- 
ion NS or BH in a close binary system. We first review the 
characteristic properties of these disks, and then discuss the 
physical approximations made in this study. 

2. 1 . Disk Properties 

As discussed in Fryer et al. (1999) and Ml 2, whether the 
WD is disrupted by the primary depends on whether mass 
transfer is stable or unstable following the onset of Roche lobe 
overflow. Stability depends on the mass ratio of the binary 
q = MwY)/M c , with higher values of q > 0.2-0.5 favoring dis- 
ruption (e.g., Paschalidis et al. 2009). Here M W d and M c are 
the mass of the WD and the central compact object (NS or 
BH), respectively. Once the WD is disrupted, conservation of 
angular momentum implies that the material will circularize 
around the NS/BH at a characteristic radius 
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is the orbital separation at Roche lobe overflow (Eggleton 
1983) and Kwd is the WD radius. Note that equation (1) ne- 
glects the self-gravity of the disk. The WD radius Rwd is 
related to its mass by (Nauenberg 1972; Fryer et al. 1999) 
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where Men — L45(/x e /2) 2 M© is the Chandrasekhar mass 
and n e is the mean molecular weight per electron. Table 1 
gives the circularization radii for models considered in this 
paper. 
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As matter circularizes around the NS/BH, a fraction of the 
kinetic energy of rotation is thermalized. This results in a 
thick torus (e.g., Fryer et al. 1999) with a scale-height Ho <~ 
Ro/2 and an average density 
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Assuming that the internal energy e; nt is dominated by non- 
degenerate particles, and that it balances 25% of the gravi- 
tational energy e grav at the circularization radius, one finds a 
characteristic torus temperature 
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is the mean molecular weight, with Y e the electron fraction, 
and {Xj,A,} the mass fraction and atomic number of ionic 
species i, respectively. 

The characteristic timescale for matter to accrete may be 
estimated by the viscous time 5 
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where a parametrizes the disk viscosity (§2.2), resulting in a 
characteristic accretion rate 
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For the thermodynamic conditions above, the opacity is 
dominated by electron scattering. The timescale for photons 
to diffuse out of the disk midplane is then 
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where K es is the electron scattering opacity. Note that since 
fdiff is much longer than the timescale for disk formation (<~ 

5 Note that ( v i sc evaluated at large radii in the disk generally underestimates 
the true accretion timescale onto the central object since the net accretion 
rate usually decreases with decreasing radii (e.g., M oc r) in most models 
of RIAFs. Equations (8) and (9) nevertheless represent useful characteristic 
values. 



forb) or viscous evolution (~ f v ; sc ), this implies that the disk 
is highly radiatively inefficient and hence radiation losses can 
be neglected. 

We limit our simulations to regions of the disk external 
to the radius R[ n ~ 10 7 cm ~ 0.01/?o, a choice made pri- 
marily for computational expediency. This cutoff is justi- 
fied because nuclear burning primarily occurs exterior to this 
radius, with only endothermic dissociation happening inside 
this point (M12). 

Neutrinos can in principle also cool the disk. The timescale 
for neutrino cooling near the disk circularization radius Rq due 
to e~/e + pair annihilation is estimated to be f„ ~ 5 x 10 s yr at 
the fiducial densities and temperatures given in equations (4) 
and (6), again much longer than the timescale for disk evolu- 
tion. However, as matter flows inwards to smaller radii where 
the temperature is higher, neutrino cooling is increasingly im- 
portant, possibly even causing the disk to transition to a geo- 
metrically thin configuration (e.g., Di Matteo et al. 2002). A 
straightforward calculation shows that a thin disk obtains in- 
side a critical radius R v given by (Chen & Beloborodov 2007; 
see also Metzger et al. 2008, their eq. [11]) 
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If the characteristic accretion rate near the outer edge of the 
disk (9) is similar to that at smaller radii in the disk (as may 
not be valid for RIAFs), then from equation (11) one infers 
that a thin disk is only possible at radii <C 10 7 cm. Although 
neutrino cooling is unlikely to significantly alter the disk ther- 
modynamics in our computational domain, we nevertheless 
include this effect in most calculations (§3.3) to compensate 
for the exclusion of radiation pressure (§2.3). 

2.2. Angular Momentum Transport 

Transport of angular momentum in a fully ionized disk is 
thought to be mediated primarily via Maxwell stresses asso- 
ciated with the MRI. If the initial magnetic field of the WD 
is strong and the torus forms with a strong poloidal field, 
then a magnetocentrifugal wind can also carry away angular 
momentum (e.g., Blandford & Payne 1982; Stone & Pringle 
2001). The interior magnetic fields of WDs are not well 
constrained observationally, but the measured exterior fields 
range from < 10 kG in most WDs, up to several hundred MG 
in a small population of highly magnetized WDs (Wickramas- 
inghe & Ferrario 2000). These fields are much less than the 
value ~ 10 12 G that would reach equipartition with the gas 
pressure in the torus, 6 suggesting that angular momentum loss 
to winds may be neglected. 

Assuming that the magnetic field is weak, the evolution of 
the torus should resemble that found by previous three dimen- 
sional MHD simulations of RIAFs onto black holes at large 
distances from the event horizon (Hawley et al. 2001; Igu- 
menshchev et al. 2003). Unfortunately, MHD simulations are 
computationally expensive. One reason is that toroidal mag- 
netic fields amplified by the MRI tend to rise buoyantly into 
low density regions above the disk midplane (e.g., Davis et al. 
2010; Shi et al. 2010). This causes the Alfven speed to be very 
large, thus limiting the dynamic range in radii or evolution 
time that can be explored due to a small Courant-Friedrichs- 
Lewy (CFL) timestep (Stone & Pringle 2001). This constraint 

6 The large-scale field may be amplified somewhat during disk formation 
via linear field winding, but this is only likely to enhance the toroidal field. 
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becomes especially severe given the need to resolve the most 
unstable mode of the MRI for a reliable saturation amplitude 
(e.g., Hawley et al. 1995). 

A less costly alternative, which we adopt, is to perform 
hydrodynamic simulations with a suitably chosen anomalous 
shear stress. There are obvious caveats to this approach. First, 
the MRI converts orbital kinetic energy directly into turbulent 
magnetic and kinetic energy, with the ensuing stresses being 
responsible for the transport of angular momentum. In con- 
trast, viscous heating converts orbital energy into heat, which 
then drives convection through entropy gradients (Hawley 
et al. 2001). This causes the hydrodynamic models to achieve 
a state of marginal convective stability (Stone et al. 1999), 
whereas in MHD the very same stability criterion must be vi- 
olated in order for the MRI to operate (Hawley et al. 2001). 
The relative spatial distribution of specific angular momen- 
tum and entropy thus differs between the hydrodynamic and 
MHD cases. In addition, the spatial and temporal distribution 
of heating is very different between MRI-driven magnetohy- 
drodynamic turbulence and shear viscosity (e.g., Hirose et al. 
2006). This could potentially alter the thermodynamic struc- 
ture of the disk and hence the spatial distribution of nuclear 
burning. 

Despite these caveats, purely hydrodynamic simulations 
capture important aspects of MHD simulations. First, the 
time-averaged mass fluxes in the disk midplane have the same 
qualitative form, with strong inflow and outflow nearly can- 
celing each other out, resulting in a net accretion rate that is 
approximately constant with radius (Stone et al. 1999; Stone 
& Pringle 2001; Hawley et al. 2001). Second, the radial scal- 
ings of time-averaged quantities (e.g., density and tempera- 
ture) found by MHD simulations can be matched by hydro- 
dynamic simulations if the functional form of the shear stress 
is suitably chosen (Stone & Pringle 2001). These scalings 
do not differ substantially between two- and three dimensions 
(Hawley et al. 2001). 

The experimental approach of this study demands the abil- 
ity to simulate accretion torii over a large dynamic range in 
spatial distances and timescales. We thus begin our study of 
NuDAFs by performing hydrodynamic simulations with an 
anomalous viscous stress. A similar approach was adopted 
recently by Schwab et al. (2012) in studying the evolution 
of accretion disks created by WD-WD mergers. In order to 
best evaluate the uncertainties introduced by this approach, 
we adopt several functional forms for the kinematic viscos- 
ity, each of which lead to different transient and quasi-steady- 
state outcomes. 

The first parameterization is that which makes the ratio of 
viscous to orbital times independent of radius, 
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where i>o is a dimensionless constant. This form was found 
by Stone et al. (1999) to yield a self-similar behavior in the 
inner regions of the disk, with p oc r" 1 ' 2 and T oc r~ l . The 
relationship between vq and the ratio of timescales is 
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We employ this prescription to compare with analytic expec- 
tations for the dynamical importance of burning (§2.3). 
We also employ a Shakura & Sunyaev (1973) parameteri- 
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where £1% is the Keplerian frequency and c 2 = jp/pis the adi- 
abatic sound speed. The functional form of v a differs from u ct 
in the additional dependence on the ratio of thermal to gravi- 
tational energies. At the circularization radius, the numerical 
value of vq in equation (13) corresponds to a ~ 0.01, with 
exact values depending on the thermal content of the disk. 

Finally, to connect with MHD calculations, we employ the 
prescription of Stone et al. (1999) 
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where p\ is a reference density, taken to be the maximum 
value in the initial condition (§3.2). This parameterization 
yields radial scalings in agreement with axisymmetric MHD 
simulations with low initial poloidal fields (Stone & Pringle 
2001). 

Following Stone et al. (1999), we include only the az- 
imuthal components of the viscous stress tensor T in our sim- 
ulations. 
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Including all components of the stress would make the fluid 
truly viscous, and suppress convection in the poloidal direc- 
tion (e.g., Igumenshchev & Abramowicz 1999). Our simu- 
lations, like previous hydrodynamic studies, mimic turbulent 
angular momentum transport via thermally-driven convec- 
tion. Suppressing this effect would suppress the very mecha- 
nism that gives plausibility to a hydrodynamic modeling. 

Given the potentially large mass of the disrupted WD rel- 
ative to that of the central NS/BH, it is also possible that the 
disk will be susceptible to gravitational instabilities, which 
could contribute to the transport of angular momentum via 
spiral density waves (e.g., Larson 1984). We neglect this pos- 
sibility here, since pressure waves may stabilize such instabil- 
ities (the disk is quasi-virial and geometrically thick). Also, 
our analysis does not include self-gravity. However, to the ex- 
tent that such waves are present and transport angular momen- 
tum outwards, some of the qualitative effects of self-gravity 
may be captured by our prescription for angular momentum 
transport. 

2.3. Microphysics 

For a given initial composition of the torus, we focus on the 
reaction activated first at large radii in the disk, since energy 
input from this reaction is the most important relative to the 
local rate of viscous heating (M12). 

We focus mostly on the reaction 12 C( 12 C,7) 24 Mg, relevant 
to disks formed from mid-mass C-0 WDs. We also calcu- 
late a few models that use the triple-alpha reaction, as would 
be appropriate for disks arising from low-mass He WDs or 
He stars. Our analysis is also applicable to a-capture reac- 
tions [e.g., 4 He( 16 O,7) 20 Ne], which are relevant for a hybrid 
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He-C-O WD disk or to the outer layers of a collapsing Wolf- 
Rayet star (Woosley & Heger 2006). We exclude the param- 
eter regime relevant to massive WDs with O-Ne composition, 
because the small circularization radii Rq (eq. [1]) and high 
disk temperatures 7^,. (eq. [6]) suggest that in these systems 
burning begins during the circularization process itself (M12). 

To better enable analytic comparisons, we adopt a generic 
power-law form for the 12 C( 12 C,7) 24 Mg reaction rate 



X f =-AXfpT' 1 



(18) 



where Xf is the mass fraction of fuel (carbon), A is a normal- 
ization constant, and /3 results from expanding the analytic 
expression for the reaction rate as a Taylor series in temper- 
ature (e.g., Kippenhahn & Weigert 1994) about the point at 
which the burning time and dynamical time (or viscous time) 
are comparable in the disk midplane. In the temperature range 
(0.6- 1.2) x 10 9 K, this procedure yields f3 = 29. In some of 
our models we employ the full functional form of each reac- 
tion (Caughlan & Fowler 1988), suitably converted to code 
units given the parametric equation of state. Details of this 
procedure are given Appendix A. 

Associated with each reaction is the specific nuclear bind- 
ing energy released, 



£nuc = QXf/m a , 



(19) 



where Q is the total energy released per reaction and m a is 
the mass of the reaction product. In several simulations we 
artificially suppress the value of Q oc e nuc from its true physi- 
cal value Qo in order to explore how the outcome of NuDAF 
evolution depends the amount of nuclear energy released. 

For simplicity, we use an ideal gas equation of state with a 
single adiabatic index 7 = 5/3. As shown in M12, gas pressure 
exceeds radiation pressure close to the disk midplane during 
phases of evolution when most of the mass accretes (i.e. on 
timescales t > f v i sc ). However, radiation pressure contributes 
an increasingly larger fraction of the total pressure at smaller 
radii in the disk. It dominates over gas pressure when the 
density is low at early times prior to when the accretion rate 
achieves a steady-state. Neglecting radiation pressure results 
in an overestimate of the temperature (and hence the rate of 
nuclear reactions), potentially producing unphysical prompt 
detonations in some of our simulations, an issue we discuss 
further in §4.2. For completeness, we also include neutrino 
cooling as parameterized by Itoh et al. ( 1996). 7 

3. NUMERICAL METHOD 

3.1. Equations Solved 

Spherical polar coordinates in axisymmetry (r,9) are 
adopted, with the origin at the center of the neutron star. We 
solve the equations of mass, momentum, energy, and chem- 
ical species conservation, with source terms due to gravity, 
shear viscosity, and nuclear reactions: 
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where p, p, \ p = v r r+\>e8, and e ml are the fluid density, pres- 
sure, poloidal velocity, and internal energy, respectively, and 
d/df = d/dt +v p ■ V is the Lagrangian differential operator. 

The specific angular momentum along the symmetry axis 
satisfies l z = rsin^v^, with the azimuthal velocity. The 
corresponding centrifugal force in the poloidal direction is 
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The system of equations (20)-(24) is closed with an ideal 
gas equation of state with adiabatic index 7 = 5/3, so that 
p = (7- \ )pei at . The nuclear abundances are constrained by 
baryon number conservation, Xf+X a +X x = 1, where X a and 
X x denote the mass fractions of ash and inert species, respec- 
tively. The nuclear energy generation rate Q nuc is given by 
one of the full reaction rates described in Appendix A, or by a 
power-law approximation using equations (18) and (24). Q coo i 
is the neutrino energy loss rate per unit mass. 

Throughout this paper, we express quantities in terms of 
the Keplerian velocity at the circularization radius = 
(GMc/Ro) 1 / 2 , and the orbital time at the same location, f or b — 
2 7 rfl 3 / 2 /(GM c ) 1 / 2 . 

3.2. Initial Conditions 

The initial condition is a constant entropy and specific an- 
gular momentum torus (e.g., Papaloizou & Pringle 1984) with 
uniform chemical composition. A realistic merger will pro- 
duce a torus with a more general radial distribution of angular 
momentum, since the disruption process is not instantaneous 
(e.g., Fryer et al. 1999). We ignore this complication in the 
interest of clarity. The implications of this assumption in light 
of our results are discussed in §5.1. 

Following Stone et al. (1999), the torus density is normal- 
ized to its maximum value p max , thereby fixing the poly tropic 
constant in terms of the adiabatic index and the torus distor- 
tion parameter d. The resulting initial density distribution is 



Pmax 



7-1 



2d 
d^l 



r 



1 ( Ro V 1 
2\rsm6J 2d 



(26) 



The distortion parameter d is a measure of the internal en- 
ergy content of the torus. The maximum ratio of internal to 
gravitational energy occurs at the point of maximum density 
(r=R , z = 0), 

^int.max 1 «~ ( r )l\ 

GM/Ro ~ 2^~d~' 

Given the angular momentum distribution j(r) of the torus, 
the value of d is in principle fully determined by the prop- 
erties of the disrupted binary and energy conservation. Un- 
fortunately, the uncertainty in j(r) and the complication of 
self-gravity (which is not included in the initial torus solution) 
impede a precise determination of d. In the absence of more 
detailed information, we adopt d = 1 .5 as a fiducial value for 
most of our models. For 7 = 5/3, this yields a maximum ratio 
of internal to gravitational energy (eq. [27]) of 10%. We also 
use d = {1.2,3} in some models, to study the sensitivity of 
our results to this choice. The integrated mass distributions so 
obtained bracket the results of Fryer et al. (1999) for WD-BH 
disks. 
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A physical temperature is obtained by assuming that the 
pressure is provided by a non-relativistic ideal gas, and using 
a mean molecular weight (eq. [7]) appropriate to the initial 
WD composition. Similarly, a physical density is obtained by 
multiplying the average disk density (eq. [4]) by the ratio of 
the maximum to average density in the initial torus. This ratio 
is a function of d and 7 only. Details of this unit conversion 
are provided in Appendix A. 

The torus is surrounded by a low-density isothermal atmo- 
sphere, 

f GM (r m \\ 
pAr) = Pat,o exp i -j — ( 1 V , (28) 

I c at r in ^ r ' ) 

where r m is the inner radial boundary of the computational do- 
main, and c at = /? a t,o / Pat,o is the isothermal sound speed of the 
atmosphere. This surrounding medium has the advantage of 
being stably-stratified. We find that a constant density atmo- 
sphere, such as that used by Stone et al. (1999), generates ex- 
cessive numerical noise near the inner boundary. The normal- 
ization constant is chosen to be much smaller than the bulk of 
the torus density, and an order of magnitude below the density 
at which we cut off source terms (3.3). At radii r ^> r m , the 
density asymptotes to p^, = j o at exp[-GM c /(c^ t ri n )]. In most 
models, we adopt p a t,o/Pmax = 10~ 6 and c 2 t = 25GM c /Rq. The 
temperature is chosen so that the resulting density scale height 
at the inner boundary is resolved with at least 5 cells in radius. 
This choice of atmospheric parameters causes a negligible ef- 
fect on the torus evolution even though the atmosphere itself is 
unbound. In a few models, we increase the background den- 
sity to p a t,o//Omax = 10~ 4 to obtain smoother evolution of the 
inner edge of the torus during the initial infall. This higher 
value artificially slows down ejected material, however, and 
is kept for comparison purposes only. 

3.3. Time-dependent evolution 

We use FLASH3.2 (Dubey et al. 2009) to evolve the system 
of equations (20)-(24) with the dimensionally-split version of 
the Piecewise Parabolic Method (PPM, Colella & Woodward 
1984). The public version of the code has been modified to 
allow for a non-uniformly spaced grid in spherical polar coor- 
dinates as described in Fernandez (2012). 

The specific angular momentum is included as an advected 
scalar subject to viscous source terms (eq. [22]). To take ad- 
vantage of the finite-volume character of the hydrodynamic 
method, the diffusion operator is recast as 

rsin0(V-T) = V-F £ , (29) 

where 

F^=rsin0(V r ^r+2V) (30) 

is the angular momentum flux vector. This diffusive flux is 
added to the advective flux obtained from the Riemann prob- 
lem solution at each cell interface. The combined flux is then 
used to update the volume-averaged value of L (see Lindner 
et al. 2010 for a similar implementation in cylindrical coordi- 
nates). The internal energy is updated by simple operator-split 
addition of the scalar viscous energy generation rate (equa- 
tion 23). Numerical stability requires the evolution time step 
to be lower than 

1 /Ar 2 r 2 A0 2 \ 
Af visc =-min , , (31) 

2 y v v J 



where the minimization is carried out over the entire compu- 
tational domain. This constraint on the time step is imposed in 
addition to the CFL condition required by the hydrodynamic 
method. The centrifugal force is included as part of the de- 
fault treatment of 'fictitious' forces that arise in curvilinear 
coordinates. 

By default, FLASH3.2 evolves the total specific energy (in- 
ternal plus kinetic). To avoid adding the additional source 
terms needed to account for the Eulerian rate of change of 
the rotational kinetic energy, we only evolve internal energy 
(equation 23). 

The nuclear energy generation rate and neutrino cooling are 
also included via simple operator-split update of the internal 
energy. We suppress nuclear burning inside shocks, and em- 
ploy a standard treatment of multi-species advection (Plewa & 
Miiller 1999). To avoid artificial thermonuclear runaways, we 
impose an additional constraint on the time-step from nuclear 
burning 

Af nuc = 0.84^. (32) 

In practice, the simulation time-step is dominated by the CFL 
condition given the atmospheric temperature and the grid res- 
olution, being usually 10 to 100 times smaller than Af nuc . 

The computational grid is logarithmically spaced in the ra- 
dial direction, covering four orders of magnitude in length. In 
most cases, it extends from r m i n = 0.01i?o to r max = 100/?o, al- 
though a few models use a smaller inner radius (Table 1). In 
the polar direction the grid covers the full range of latitudes 
([0°,180°]), with cells having a constant ratio of sizes (e.g., 
Stone & Norman 1992). We make this spacing symmetric 
relative to the equator, with the coarsest cells next to the polar 
axis. Our standard resolution is N r = 64 cells per decade in ra- 
dius. The angular spacing is chosen so that Ad ~ Ar/r ~ 2° 
at the equator, and A9 p = 4° or 6° next to the polar axis 8 . 
We evolve one model at higher resolution, with N r = 128 cells 
per decade in radius, AO ~ Ar/r ~ 1° at the equator, and 
A9 P = 2°. 

The boundary conditions at the polar axis are reflecting in 
all variables. At the inner radial boundary, we impose a zero- 
gradient boundary condition for the mass flux, subtracting the 
contribution from the isothermal atmosphere. For a ghost cell 
with position (r, 8), the density, pressure, and radial velocity 
are set to 

p(r,9) = p(r u 9)-p at (r l ) + p, l (r) (33) 
p(r,9) = p(r u 9)-p, t (r l )+p, t (r) (34) 

v r (r,6)=(r l /r) 2 v r (r u 9) (35) 

if v r (r\ , 9) < and p(r\ , 9) > p a t('"i X where r\ is the radial co- 
ordinate of the center of the first active cell outside the inner 
radial boundary. Otherwise the radial velocity in the ghost 
cells is set to zero. Because fluid elements associated with 
the disk normally have densities much larger than that of the 
atmosphere, this is in effect a standard outflow boundary con- 
dition. The additional terms ensure that the isothermal atmo- 
sphere remains quiescent when no material is accreted. The 
ghost cells for the internal energy are set in consistency to 
those of the density and pressure. All other variables are given 
a zero radial gradient in the ghost cells. This prescription is 
repeated at the outer radial boundary, with the only modifica- 

8 We find that the funnel that develops around the rotation axis can cause 
numerical problems if the angular resolution is too fine. 
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TABLE 1 
Models Evolved 



Model 


M\VD 

(Mq) 


XH c /Xc/Xo' d 




Reaction 


e/Go d 


Visc. e 




/V r s 




>"m/Ro 


Rq 

(cm) 


Poo' 
(Pmax ) 


Po J 


v Cool? 


COqOOO 


0.6 


0/0.5/0.5 


1.5 


power-law 


0.00 


ct 


0.0016 


64 


4 


lO" 2 


10 9 - 3 


IO" 57 


10~ 3 


Yes 












U. 1U 






























0.25 




















COq050 










0.50 




















COq075 










0.75 




















COqlOO 










1.00 




















COq000_nc 


0.6 


0/0.5/0.5 


1.5 


power-law 


0.00 


ct 


0.0016 


64 


4 


10 - 


1 rt*3 T 
10 


i A-7 7 

10 


10 3 


No 


COq010_nc 










0.10 




















V,UL|Ui: J I1C 






























COq050_nc 










0.50 




















COql00_nc 










1.00 




















CUqU50_HK 


0.6 


0/0.5/0.5 


1 .5 


power-law 


0.50 


ct 


0.0016 


128 


2 


10 


10 


10 


10 


Yes 


Cuq050_at 
















64 


4 






i a-7 7 
10 






COq050_p 


























10^ 




COq050 dL 






1.2 




















io- 3 




COq050 dH 






3.0 












6 












COq050_i/H 






1.5 








0.016 
















COq050_spb 












spb 


0.008 




4 












COq050_a 












a 


0.0144 
















CO fl 


0.6 


0/0.5/0.5 


1.2 


i2 C+ 12 C 


1.00 


ct 


0.016 


64 


4 


io- 2 


to 9 - 3 


io-" 


io- 3 


Yes 


CO_f2 






1.5 








0.0016 
















HE 00 


0.3 


1/0/0 


1.5 


3a 


0.00 


ct 


0.016 


64 


4 


io- 25 


10 96 


10 -11.8 


io- 3 


Yes 


HE_fl 










1.00 





















a Initial mass fractions of helium, carbon, and oxygen, respectively. 

b Distortion parameter controlling the initial thermal content of the torus (see eq. [27]). 

c Nuclear reaction rate employed (Appendix A). 

d Ratio of the nuclear energy released per reaction Q to the true physical value Qq. 
e Functional form of the kinematic viscosity (eqs. [12]-[15]). 

f Magnitude of viscosity i>o or a, as appropriate. The coefficient a is related to uq through equation (27). 

g Number of cells per decade in radius. 

h Angular resolution at the polar axis 

1 Asymptotic atmospheric density (eq. [28]). 

J Cutoff density for source terms (eq. [36]). 



tion being a reversal of the required sign of the radial velocity 
for outflow. 

To prevent excessive heating in regions of low density, we 
multiply all the energy and viscous source terms by a cutoff 
function f(p). Otherwise the CFL timestep can become pro- 
hibitively short in regions that are unimportant for the disk 
evolution. The functional form of the cutoff is 

where po is a fiducial density, which we choose to be 
Po/ Pmax = 10~ 3 . In the initial torus configuration, 99.99% 
of the mass is included inside this contour. We have veri- 
fied that the evolution of the relevant dynamic quantities is 
not very sensitive to the value of this cutoff below a certain 
value (§3.4). Although a low-density suppression of burning 
may seem ad hoc, it can be motivated physically. First, at 
low densities our assumed ideal gas EOS overestimates the 
true temperature (and hence burning rate), since we have ne- 
glected the effects of radiation pressure, which dominates the 
total pressure in low density regions. Also, in low density re- 
gions the photon optical depth may become sufficiently small 
that radiative cooling (also neglected in our simulation) may 
become relevant on the timescales of interest. 

Numerical tests of the angular momentum evolution are 



presented in Appendix B. 



3.4. Models Evolved 

We evolve a suite of models that systematically sample the 
parameter space of disk properties relevant to WD-NS merg- 
ers. The complete set is shown in Table 1. 

Models are grouped according to the parameter being var- 
ied. Most cases focus on a C-O WD of mass Mwd = 0.6Mq 
accreting onto a NS of mass M c = 1.4M Q . The composi- 
tion is 50% carbon and 50% oxygen. The first two subsets 
of models focus on the effect of varying the energy released 
per nuclear reaction Q, expressed as a fraction of its true 
physical value Qq, using the power-law approximation to the 
12 C( 12 C,7) 24 Mg reaction (eq. [18]). One sequence includes 
neutrino cooling and the other does not (models ending in 
_nc). 

As we discuss below, most of our models with Q > 0.5(2o 
undergo large-scale detonations, while those with Q < 0.5go 
instead undergo quiescent burning. Since one of our goals 
is to understand what conditions are necessary for detona- 
tion, we adopt the marginal case of the sequence with cool- 
ing (COq050) for further parameter variation. The third sub- 
set of models thus explores the effect of higher resolution 
(COq050_HR), torus distortion parameter (COq050_dL and 
COq050_dH), strength (COq050_!/H) and functional form 
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FIG. 1 . — Snapshots of model COq050_HR, illustrating the main evolutionary stages of the accretion of a tidally disrupted 0.6Mq C-0 WD onto a 1.4Mq 
NS in a case that leads to a detonation. The mass fractions of fuel (carbon) and ash (magnesium) are shown in blue and yellow, respectively, with green regions 
indicating mixed material. The shading is proportional to the density gradient. From left to right, panels correspond to (a) the initial condition; (b) establishment 
of a steady burning front at r nu c ~ 0.1i?o (e<3- [39]); (c) initiation of detonation; and (d) shock expansion and partial unbinding of the torus. Time is shown in units 
of the orbital time at r = Rq. Note that the spatial scale is enlarged in the final two panels. An animated version of this figure is available in the online version of 
the article. 



(COq050_spb and COq050_a) of the shear stress, magnitude 
of the atmospheric density (COq050_at), and value of the cut- 
off density po (COq050_ / o0) on the marginal model. 

The fourth subset of models employs the full reaction rates 
to assess the likelihood of detonation in disks formed from 
C-O and He WDs. Model CO_fl employs the same param- 
eters as our fiducial C-O WD marginal model, but with the 
full strength of the energy release (Q = Qo). Model CO_f2 de- 
creases d and increases vq to increase the chance of explosion. 
Models He_000 and He_f 1 explore the evolution of pure He 
disks (^He = 1 ) under the influence of the triple-alpha reaction. 
Parameters correspond to a Mwd = 0.3M Q WD accreting onto 
a NS of mass M c = 1.4M©. The viscosity is chosen large, to 
maximize the likelihood of a detonation. The inner boundary 
of the domain is chosen smaller than in the C-O case, because 
the ratio of burning to circularization radii is expected to be 
smaller as well (M12). 

4. RESULTS 

We begin the presentation of our results with a basic 
overview (§4.1) of the impact of nuclear burning on the evo- 
lution of the accretion disk, focusing on a model which deto- 
nates. Then, in separate sections, we elaborate on exploding 
(§4.2) and quiescent cases (§4.3). The results of our simula- 
tions are summarized in Table 3. 

4.1. Overview 

In the absence of nuclear burning, the disk evolves in a way 
similar to that in the axisymmetric simulations of RIAFs by 
Stone et al. (1999). The initial torus develops a radial velocity 
through the action of viscous stresses. Material near the inner 
edge of the disk reaches the inner boundary within the first 
few orbital periods as measured at r = Ro- As more of the disk 
spreads inward, viscous heating drives convective turbulence, 
resulting in inward and outward mass fluxes of comparable 



magnitude. Material is ejected in a wide funnel around the 
rotation axis, but most of the ejecta has negative energy and 
eventually falls back to the disk. The accretion rate peaks on a 
timescale which is comparable to a viscous time at r = Ro, be- 
fore gradually decreasing with time thereafter. The properties 
of torii without burning are discussed in relation to quiescent 
NuDAFs in §4.3. 

The dynamical importance of nuclear energy generation at 
a given radius can be quantified by the ratio of the specific 
nuclear energy released per reaction e nuc (eq. [19]) to the lo- 
cal gravitational binding energy. Given the steep temperature 
dependence of most reactions, fuel depletion and energy de- 
position are generally concentrated in a narrow surface that 
crosses the midplane at a characteristic radius r nuc from the 
central object. If burning occurs subsonically, the magnitude 
of the density jump across the burning region is given by the 
ratio of e nuc to the specific enthalpy of the fluid w (Landau & 
Lifshitz 1987; see also Appendix C). The latter is related to 
the specific gravitational binding energy by a factor oc (H / f) 2 
(e.g., M12). Hence, one can use the ratio of e nuc to w evalu- 
ated at r = r m , r , 



w(r nuc ) 
2d 



d- 1 GM C /R 



' nuc 

~Ro~ 



(37) 



(38) 



to parameterize the importance of nuclear reactions on the 
disk evolution. Note that equation (38) is specific to our ini- 
tial conditions (§3.2); it assumes that the temperature has a 
power-law dependence with radius, T oc r~ m , and makes use 
of equation (27) for normalization. The values of 'J quoted in 
our results below are calculated using an angle-averaged value 
of w(r nuc ) obtained directly from our simulations. Equation 
(38) is a reference value only. 
The location of r nuc can be estimated by equating the dy- 
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TABLE 2 

Dynamical Importance of Nuclear Reactions in RIAFs. a 



Reaction 


2o 
(MeV) 




* nuc 

(10 8 cm) 


v]> d 


12 C( 12 C, 7 ) 24 Mg 


13.93 


0.0/0.5/0.5 


1.9 


1.68 


12 C(q,7) 16 


7.16 


0.2/0.5/0.3 


2.5 


1.74 


16 O(«,7) 20 Ne 


4.73 


0.2/0.3/0.5 


5.1 


1.86 


16 0( I6 0,7) 32 S 


16.54 


0.0/0.5/0.5 


1.0 


0.78 



r l l 2 (M oc r); viscosity parameter a = 0.01; H/r = 0.3, and M c = 1.4Mq. 
b Initial mass fractions 

c Burning radius at which half the initial fuel is consumed. Its value is calcu- 
lated numerically using the full functional form of the reaction rate. 

Ratio of £nuc to the enthalpy w = (7 - \T l (H/r) 2 GM c /r mlc of the flow 
(eq. [37]). 

namical time r/Vko to the nuclear burning time e mt /Q nuc . If 
the timescale for viscous heating is long compared to that 
for vertical hydrodrostatic balance (the sound crossing time), 
then the disk temperature T approaches its virialized value 
r v ; r oc r~ l (eq. [6]). Otherwise it has a more general depen- 
dence T oc r~ m . If the density scales as a power-law with ra- 
dius p oc r", then using the power-law approximation to the 
12 C( 12 C,7) 24 Mg reaction (eq. [18]), and the power-law form 
for the viscosity prescription (eq. [12]), one finds that 



Ro 



^nuc 
eint 



AXfp T£ 



l/(m[/3-l]+n-3/2) 



(39) 



where the zero subscripts on thermodynamic quantities de- 
note their values at r = Rq. For C-O torii with Mwd = 0.6M Q , 
M c = l.4M@, and d = 1.5, we find r nuc /Ro ~ 0.1, for m = 1, 
and n = 2. Equation (39) is sensitive to the radial dependences 
of the temperature and density, as well as to the temperature 
normalization at the circularization radius. The position of 
r nuc will thus change in time as the disk approaches a quasi 
steady-state configuration. 

Assuming a carbon mass fraction Xq = 0.5 in e nuc , a value 
of r nuc ~ 0.1 yields "J/ ~ 1.8(<2/<2o) for the power-law burn- 
ing rate and a central mass M c = 1.4M Q . Using the full func- 
tional form for the reaction rate requires solving for r nuc /^o 
numerically. Table 2 shows the result of such a calcula- 
tion for 12 C( 12 C,7) 24 Mg and other relevant reactions (e.g., a- 
captures). Note that ^ is similar in all cases except for oxygen 
burning, which is lower by a factor ~ 2. 

The qualitatively new behavior introduced by nuclear burn- 
ing is illustrated in Figure 1, which shows several snapshots 
in the evolution of the C-O torus in our fiducial high resolu- 
tion model COq050_HR with Q = 0.5Q (* ~ 0.6 measured 
at the moment of final detonation, see Table 3). Initially, the 
disk is composed entirely of fuel (carbon). During the first 
orbital period of evolution, the inner edge of the disk moves 
inwards and increases in temperature, before igniting at some 
radius initially smaller than the fiducial burning radius r nuc . 
A narrow burning front is established which moves outwards 
and then settles at r ~ r mc on a (local) dynamical time. Tur- 
bulence is generated outside of the burning front, triggering a 
thermonuclear runaway and detonation at the time t = 1 .67f 01 -b. 

After a time t = 1.9f or b, the detonation has successfully 
propagated out through the entire torus. At this point approx- 
imately 99% of the initial WD mass (~ 0.6M Q ) has become 
unbound. The net energy in this material is E e ; ~ 3 x 10 50 



U, l 




r/R 

FIG. 2. — Panels (a) and (b): Angle-averaged density p and tempera- 
ture T as a function of radius for the model COq()50_HR (Fig. 1), re- 
spectively. Averages are density-weighted and taken within 45 degrees of 
the equator. From blue to red, solid curves correspond to times r/( or b = 
{0,0.5, 1, 1.2, 1.673, 1.68, 1.69, 1.7}. The horizontal dashed gray line shows 
the cutoff density po below which nuclear burning is artifically suppressed 
(eq. [36]). Panel (c): Various timescales in the disk as a function of ra- 
dius at the time f /f rb = 1.673 just prior to detonation (compare with Fig. 3). 
Timescales f, shown include: dynamical fj yn = r/v^if) (grey); sound cross- 
ing 'sound = r/c, (black); viscous heating fv isc -heat = eint/Gvisc (blue); nuclear 
burning fburn-heat = <?int/2nuc (red); and neutrino cooling kern-heat = eint/Gcool 
(green). 



ergs, and has a mass-weighted radial velocity of v e 



6000 



km s" 1 . A fraction 5% of this mass becomes bound and re- 
mains as a remnant disk. In §5 we discuss the possible obser- 
vational signatures of this supernova-like explosion. 

Our calculations show that detonations can occur in 
NuDAFs ranging from microexplosions (e.g., Lisewski et al. 
2000) which remain localized, to complete consumption of 
the initial disk. A large-scale detonation capable of unbinding 
the disk, and occurring within the first few orbital periods of 
evolution, appears to be a robust outcome if "J/ is above a criti- 
cal value "fait ~ 1 . A primary cause of such detonations is that 
the shape of the initially stationary burning front is distorted 
due to the effects of the Rayleigh-Taylor (RT) instability (e.g., 
Bell et al. 2004) and convective turbulence in the surrounding 
fluid. In Appendix C we argue that the mixing of fuel and ash 
due the convective motions present in RIAFs is indeed suf- 
ficient for a detonation to occur via the Zel'dovich gradient 
mechanism (Zel'dovich et al. 1970) if 'J 7 has the appropriate 
magnitude. 

If \& <C ^cnt. then large scale detonations do not develop 
and accretion remains relatively steady across the burning 
front. During this 'quiescent' mode of NuDAF evolution, 
some fraction of the ash behind the burning front is accreted, 
while the rest is ejected into an outflow along the polar axis. 
Even when a detonation does not occur, NuDAFs produce 
a more powerful outflow than in otherwise identical RIAFs 
without nuclear burning. The fraction of the material which 
is gravitationally unbound versus that being accreted is found 
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FIG. 3. — The effect of turbulence near the burning front on the onset of detonation in the model COq050_HR. Each panel shows snapshots of the temperature, 
with the time in units of the orbital time at r = Rq displayed in the upper left corner. White contours correspond to fuel mass fractions Xf = {0.1,0.4}. From 
left to right, the panels show: quasi-steady-state burning front (f/f 01 -b = 1.600), formation of an eddy by turbulence (f / f or b = 1.650), hot spot prior to burning 
( f /'orb = 1-673), hot spot after burning (f/f or b = 1.764), and propagation of detonation (f/; or b = 1-675 and 1.679). An animated version of this figure is available 
in the online version of the article. 
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FIG. 4. — Radial profile of quantities angle-averaged over the hot spot that 
triggers a detonation at time t = 1.673f or b in the model COq()50_HR (Fig- 
ure 3). Shown are the mass fraction of fuel (green), ash (blue), the sound 
speed (black), and the inverse of the induction time gradient (red). Mixing of 
ash upstream of the front is responsible for the increased temperature. This 
heating results in the condition dr/d?j n( | ~ c, for spontaneous detonation (Ap- 
pendix C) being approximately satisfied across the hot spot. The induction 
time is defined here as (j nI j = e mt /[f3 Qnuc], where /3 = 29 is the temperature 
exponent of the reaction rate (eq. [18]). 

to increase with \I> (§4.3, Fig. 7). 

4.2. Detonations 

We begin by elaborating on the evolution of torii which un- 
dergo detonation, basing our discussion on the high-resolution 
model COq050_HR (Figure 1). The angle-averaged thermo- 
dynamic quantities in the disk as a function of radius are 
shown in Figure 2 at different times. Also shown are vari- 
ous timescales in the disk as a function of radius, calculated 
at a time t = 1 .673f rb, corresponding to the onset of the final 
detonation. 

Because the sound crossing time is everywhere shorter than 
the viscous heating time, pressure equilibrium rapidly estab- 
lishes a temperature profile that is nearly virial (T oc r" 1 ) in 
regions of the disk where neutrino cooling is negligible. This 
profile becomes shallower than r~ l interior to the burning ra- 
dius, because neutrino cooling becomes important and offsets 
viscous heating. The density also obeys a power-law radial 



profile (for p > po; eq. [36]), with a slope that decreases with 
time as matter flows inwards. Such a profile steeping is ex- 
pected since (at least in regular RIAFs) one expects p oc r" 1 / 2 
once steady-state accretion is reached on a timescale t ~ f v j sc . 

The inner edge of the torus ignites first, starting interior to 
the burning radius r nuc predicted by equation (39). Burning 
requires a significantly higher temperature (and hence smaller 
radii) than at later times because initially the density is lower 
than po (eq. [36]), below which nuclear reactions are artifi- 
cially suppressed. If neutrino cooling is not included in our 
calculations, then a global detonation is triggered at these 
very early times if the value of ^ is sufficiently high (model 
COql00_nc). Given that neutrino cooling is physically moti- 
vated and that detonations via the Zel'dovich mechanism may 
be suppressed if the effects of radiation pressure were prop- 
erly included (Appendix C), we believe that such prompt det- 
onations are unphysical, at least within the parameter regime 
of this study. 9 Verifying this will require simulations using a 
more detailed model of the initial structure of the disk and a 
more physical EOS. 

When neutrino cooling is included or if 'J 7 is lower, then a 
prompt detonation is avoided and a steady-state burning front 
develops. The front settles at the expected radius r r nuc 
within a few local dynamical times. A detonation can still be 
triggered at later times as the result of turbulence generated 
by fluid instabilities if ^ is sufficiently high. The RT insta- 
bility is somewhat suppressed by neutrino cooling when the 
value of ^ is sufficiently large. At low 'J, or in the absence 
of cooling, the burning front is noticeably distorted. Turbu- 
lence should be a ubiquitous feature of a more realistic MHD 
model since it can be generated directly by the MRI, rather 
than originating second-hand from convective instabilities. 

Figure 3 shows a few snapshots of the region localized 
around the burning front in model COq050_HR, illustrating 
how turbulence triggers a detonation. Between 1 .65 and 1 .673 
orbits, a turbulent eddy near the disk midplane mixes some of 
the hot downstream ash with the 'cold' upstream fuel, creat- 
ing a hot spot just upstream of the burning front. If mixing 

9 During the tidal disruption of a more massive O-Ne WD, nuclear burning 
may begin during the circularization process itself. A dynamical detonation 
appears more likely to occur in this case. 
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FIG. 5. — Radius of the outer edge of the burning front at the disk equator 
rf r as a function of time, shown for several models that undergo explosions. 
The radius is calculated by searching for the outermost point where the ash 
mass fraction is equal to 0.25. Panel (a) shows models with high atmospheric 
density, while panels (b) and (c) show models that vary different parame- 
ters around the marginal model with low atmospheric density and neutrino 
cooling COq050_at. Note that the sawtooth shape of some curves signals 
multiple failed detonations in cases where * is close to the critical value 
\Pcrit — 1 required for a successful large-scale detonation. The burning radius 
is r nuc ~ 0. lRo for the models shown (Table 3). The dashed curve in model 
COq()50_spb in panel (c) denotes the outer edge of the torus; in this model 
the leading shock decouples from the burning front and stalls at r ~ 50Rq. 

of ash and fuel is such that that the inverse gradient in the 
induction timescale |Vfi n d| _1 across the eddy exceeds the lo- 
cal sound speed, then a localized burning front can transition 
into a detonation via the so-called Zel'dovich grandient mech- 
anism (Zel'dovich et al. 1970; see Appendix C). We believe 
this provides a reasonable explanation for the delayed detona- 
tions observed in our simulations. 

Figure 4 shows the mass fractions of fuel and ash across the 
hot spot, confirming that mixing is indeed involved in raising 
the fuel temperature. Also shown are the sound speed and the 
inverse of the induction time gradient across the eddy. The 
latter peaks at the hot spot at a value within a factor of <~ 2 
of the sound speed, roughly consistent with the Zel'dovich 
threshold for a spontaneous detonation. 

We now discuss the fate of detonations when parameters of 
the system are varied. The discussion is centered on varia- 
tions in ^> (eq. [37]), as it directly probes the dynamical effect 
of nuclear reactions on the disk. For each model, the value 
of this parameter shown in Table 3 is measured by first find- 
ing r nuc just before explosion, and then computing the ratio 
of £nuc to the angle-averaged and mass-weigthed enthalpy at 
f = r nuc- For non-exploding models, a time average value is 
taken. From equations (38) and (39) one can see that the 
(time-dependent) radial dependence of the temperature and 
density on radius have the most influence on the time varia- 
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FIG. 6. — Ratio of the mass in the torus M tor to the initial WD mass as a 
function of time, for models that explode (Figure 5). The torus is defined as 
including all material with density p > 10~ 3 pmax. Non-exploding torii are not 
shown, as the decrease in mass is negligible on the same vertical scale. After 
explosion, the ash torus resumes accreting onto the central compact object, 
but at a reduced rate compared to disks which do not detonate (see Figure 7). 



tion of ^ in a given model. Nonetheless, the values obtained 
are generally close to what can be inferred from the initial 
properties of the disk. 

Depending on the magnitude of \I> relative to ^> crlt , the hot 
spots generated by turbulence can produce three different out- 
comes. For < "Jcrit, detonations remain localized and die 
out after propagating a distance comparable to the local ra- 
dius. Multiple such microexplosions occur as the density of 
the inner accretion flow increases with time (Fig. 3). Although 
some of these temporarily enhance the thickness of the disk, 
accretion continues relatively steadily and the burning front 
remains stationary at r « r nuc . 

When is in the vicinity of ^cnt, hot spots can trigger large 
scale detonations, but many of them die out as they propa- 
gate outwards through the disk. Because burning occurs at 
a radius r nuc interior to where the density of the torus peaks 
(r ~ Ro), such detonations are required to propagate against 
the density gradient. The weakening of the shock leads to a 
decrease in the resulting overpressure, decreasing the burn- 
ing time behind the shock and hence the energy released over 
the region where the phase velocity of burning is supersonic. 
For slow enough burning, the shock degenerates into a pres- 
sure wave and damps (Niemeyer & Woosley 1997). Models 
with "J ~ \J> crit are thus marginal in that they experience multi- 
ple failed detonations before a final explosion. Success is en- 
abled in some cases by the changing radial density gradient, 
which decreases as the inner torus viscously evolves (Fig. 2). 
Since the probability of a successful detonation increases with 
time, one must follow the evolution for a sufficiently long time 
before establishing whether a marginal model will ultimately 
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a Radius of the steady-state burning front. For delayed explosions, measured just before the last detonation sets in; time-averaged 
value for quiescent models. Not measured in prompt explosions. 

b Ratio of £nuc to the enthalpy (angle-averaged and mass-weighted) at the burning radius (eq. [37]), right before explosion or 
time-averaged if non-exploding. 

Time of the final detonation (exploding models only). 

d Expansion velocity at the time when shock reaches the domain boundary (100i?o), in units of (GM c /Rq) 1 / 2 (exploding models 
with low poo only). 

c Mass of the remnant torus (only applicable to exploding models). 

f Time range for average properties (quiescent models) or remnant mass (exploding models) calculation. 
g Time-averaged net accretion rate through the inner boundary. 
h Time-averaged net mass-loss rate at r = ARq. 

1 Ratio of net time-averaged mass-loss in nuclear ash to total ejecta, evaluated at r = 4i?o- 
J Time-averaged mass-loss rate at r = 4Rq, including unbound material only. 

k Outcomes: quiescent accretion (q), delayed explosion (d), multiple detonations before explosion (m), prompt explosion (p), 
ejecta stalls (s). 



succeed or fail to explode. 

For large VP ^> 'i'crit, the energy released by even the first 
detonation is sufficient to pass the point of maximum den- 
sity. Our parametric sequence of simulations indicate that this 
critical value "fa-it 1 appears to be relatively insensitive to 
whether neutrino cooling is included or not. The final column 
of Table 3 summarizes the outcome of each simulation. 

Figure 5 shows the time evolution of the burning front ra- 
dius rf r at the disk midplane in several models which ulti- 
mately result in explosions. When ^> is large (such as in 
COqlOO), then the burning front remains near r nuc until the 
first detonation occurs. The sawtooth shape of the burning 
front evolution for 'marginal' models (such as COq050) sig- 
nals multiple failed detonations before the final successful 
explosion. The top and middle panels of Figure 5 compare 
rf r (f) in models with and without neutrino cooling. Chang- 
ing the resolution does not alter the qualitative outcome as 
inferred from comparing models COq050 and COq050_HR. 
The higher resolution model explodes earlier while still un- 
dergoing multiple detonations. 

Variations around the marginal model with low atmospheric 



density (COq050_at) show that a lower thermal content and 
faster accretion are more conducive to detonations. Models 
COq050_dL and COq050_i/H both explode earlier than the 
marginal model. A lower cutoff density po also causes the 
disk to explode earlier. Note however that the value of is 
remarkably close to 0.6 for all the exploding models that use 
the power-law burning rate. In contrast, model COq050_dH 
does not undergo a global detonation. Counterintuitively, this 
model has a higher temperature normalization, but the deto- 
nation fails to sweep through the whole disk. We surmise that 
the larger extent of the disk plays a role in this failure. 

The functional form of the viscous stress does not appear 
to greatly influence the outcome, so long as the strength is 
comparable. Models COq050_a and COq050_spb probe the 
effects of using the viscosity prescriptions in equations (14) 
and (15). The amplitude is chosen so that the inner edge of 
the torus reaches the inner boundary at around the same time 
as in model COq050_at. Whereas model COq050_a under- 
goes a vigorous explosion (so much so that the model crashes 
after 3 orbits), model COq050_spb instead develops a prompt 
explosion that fizzles because the leading shock and the burn- 
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FIG. 7. — Mass accretion rate M m at the inner boundary (r m = O.OIRo) as a 
function of time, shown for models with low atmospheric density. All models 
shown (except for COq050_at) neglect neutrino cooling. Note that as Q [\]>] 
is increased from to Qq [1.8], the peak accretion rate decreases by over two 
orders of magnitude. The decrease in M m arises either as the result of a large- 
scale detonation (when ^ > 'J'cnt) or due to enhanced convection and polar 
outflows (when VE" < ty ait ; see Fig. 8). 

ing front decouple from one another. 

Models CO_f 1 and CO_f2 employ the full functional form 
of the 12 C( 12 C,7) 24 Mg reaction, and assume the full energy 
release (Q = Qq). They differ in their thermal content and 
strength of the viscosity. The model with higher thermal con- 
tent and lower viscosity does not explode (CO_fl), whereas 
the other (CO_f2) does undergo a large scale detonation. Be- 
cause a reasonable variation in the model parameters results 
in a qualitative change in the predicted outcome, this implies 
that the true evolution is sensitive to the details of the accre- 
tion physics. Determining whether disk detonation is indeed 
a robust outcome of C-O WD-NS mergers in Nature will thus 
require simulations which include both a more physical EOS 
and account for other details, such as the true MHD nature of 
the disk turbulence. 

Finally, we find that the He disk model with the triple-alpha 
reaction (HE_f 1) fails to detonate, in part because the burning 
begins at a radius where "J is small. Also, the weaker temper- 
ature dependence of the reaction rate leads to a burning front 
that is more spread out in radius than in the C-O WD case, 
leading to a more spatially distributed energy release. This 
situation may be less conducive to detonations generated by 
mixing across a well-defined burning front. 

Table 3 gives the time of the final detonation for each ex- 
ploding model and the shock velocity once it reaches the outer 
simulation boundary at r = 100/?o- In models with low atmo- 
spheric density {p^ = 10~ 7 7 p max ) the velocity asymptotes to a 
nearly constant value of order of the Keplerian speed near the 
radius of the original torus Vko = {GM c /Rq) x / 2 (<~ 3000 km s" 1 
in our fiducial C-O WD models). By contrast, in models with 
high background density (p^ = 10~ 5 7 /9 max ) the velocity de- 
creases and the radius of the burning front comes to an effec- 
tive halt. This outcome is an artifact of the high inertia of the 
atmosphere, as is clear by comparing the evolution of mod- 
els COq050 and COq050_at, which differ only in the value of 
Poo. In order to accurately predict the asymptotic properties 
of outflows from the disk, and hence their resulting observa- 
tional signatures, one must use a sufficiently low atmospheric 
density for which the shock motion is converged. We have 
not undertaken such a convergence study here, since our pri- 
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FIG. 8. — Top: Comparison of time-averaged mass fluxes as a function of 
radius in non-exploding models COq000_nc (no burning, black curves) and 
COq025_nc (Q/Qo = 0.25, cyan curves). Outflow, inflow, and net fluxes are 
shown as dot-dashed, dashed, and solid curves, respectively. The vertical 
dotted line corresponds to the time-averaged burning radius r mi: . Bottom: 
Root-mean-square Mach number as a function of radius in the disk midplane 
for the same set of models. The energy deposition from nuclear burning 
causes a localized increase in the turbulent kinetic energy, decreasing the net 
mass flux through the disk inside fnuc 

(compare with Figure 7). 

mary goal is understanding the basic dynamics of the disk on 
smaller scales. Also note that the density suppression of the 
burning rate (eq. [36]) ensures that no burning takes place af- 
ter the detonation has swept the initial torus material. 

Disks that explode leave behind a remnant torus composed 
primarily of ash. Figure 6 shows the time evolution of the 
torus mass M tor , defined as that enclosed within the density 
contour p > 10~ 3 p max , in several disk models. The remnant 
disk after explosion contains between 1% and 10% of the ini- 
tial WD mass, depending on the value of and on whether 
neutrino cooling is included. The existence of such a disk im- 
plies that the accretion rate onto the central compact object 
does not immediately fall to zero in disks that explode (see 
also Figure 7). 

4.3. Quiescent Disks and Explosion Remnants 

Although small-scale runaway burning and localized mi- 
croexplosions may occur for vT/ < ^ait, the disk as a whole 
never detonates. The properties of disks with nuclear heat- 
ing are nevertheless qualitatively different from those without 
burning. In this section we discuss these NuDAFs which un- 
dergo quiescent accretion. 

Figure 7 shows the time evolution of the net mass accretion 
rate at the inner boundary M- m for the sequence of C-O torii 
without neutrino cooling, for various values ^ . This sequence 
was chosen since its low atmospheric density p^ allows us 
to better study the properties of disk outflows. The asymp- 
totic accretion rate mono tonic ally decreases with increasing 
with M m dropping by more than two orders of magnitude 
as 4" increases from to ~ 1 . For large values of \I> > Vtcrit, 
this decrease in the accretion rate is the result of the mass 
loss caused by the disk detonation, which leaves only a small 
remnant torus (Fig. 6). However, even for disks that do not 
explode, the mass accretion rate is substantially reduced by 
the effects of nuclear heating (e.g., model COq025_nc). 
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contributed by nuclear heating. The light blue contour in the energy panel corresponds to pj p raax = 10 -3 , marking the approximate surface of the torus. The dashed 
white contour shows the region where the mass fraction of atmospheric material is more than 0.1% by mass (the atmosphere is very unbound by construction; 
we set its energy to zero in this plot to maximize contrast between regions dominated by torus material). The solid white contour on the heating ratio panel 
corresponds to Qvisc = 10 v|n/^0> while the dashed contour denotes the burning front (X a = 0.25). Note that both heating terms are suppressed for p < 10~ 3 p max 
(eq. [36]), hence there is negligible heating inside the polar funnel. Also note the different spatial scale of the rightmost panel. 
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gle, for three quiescent models that differ in the strength of nuclear burning. 
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nuclear ash. 



The decrease in the net accretion rate in quiescent NuDAFs 
results from more vigorous convection, driven by the en- 



hanced heating from nuclear reactions. This is illustrated in 
Figure 8, which shows the time- and angle-averaged mass 
fluxes as a function of radius for models with and without 
nuclear burning. A localized enhancement in the turbulent 
kinetic energy (as measured by the root-mean-square Mach 
number) around the time-averaged burning radius coincides 
with a sharp drop in the mass fluxes in model COq25_nc. 

For comparison, Figure 7 also includes a low-poo model 
with neutrino cooling included (COq050_at). The net accre- 
tion rate for this model is a smoother function of time than 
the otherwise identical model without cooling (COq050_nc), 
due to the near cancellation of viscous heating by neutrino 
cooling (c.f. Figure 2, bottom panel), which suppresses con- 
vection. This is also responsible for the larger magnitude of 
the net mass flux. 

In order to better understand the properties of disk outflows 
from quiescent NuDAFs, we perform time-averages of the 
flow over several orbits once they have reached a quasi-steady 
state (c.f. Stone et al. 1999). Figure 9 shows two-dimensional 
maps of density, total energy, ash mass fraction, and heating 
strength for model COq025_nc, averaged between orbits 1 1 
and 16. Material is unbound when its total energy 
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is positive. Figure 9 shows that disk material satisfying this 
condition is confined to within a narrow funnel ^10° from 
the polar axis. This unbound flow is composed primarily of 
nuclear ash and inert species. The distribution of nuclear ash 
is wider than this funnel, however. Figure 9 also shows that 
a significant fraction of the outflowing ash is bound, circulat- 
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FIG. 1 1. — Comparison of time-averaged properties between a model with nuclear burning (COq025_nc) and one without (COqO()0_nc). The average is taken 
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ing above the disk and returning to the midplane at distances 
larger than the outer radius of the disk. The angular distribu- 
tion of the total energy, radial velocity, and ash for the quies- 
cent models with no neutrino cooling are shown in Figure 10. 
Note that the expansion velocity along the funnel can reach 
the Keplerian velocity Vko at the circularization radius, which 
for the fiducial C-O WD-NS merger is ~ 3000 km s _1 . 

The rightmost panel in Figure 9 shows the fraction of the 
total heating in the disk contributed by nuclear burning. Nu- 
clear heating is important relative to viscous heating in a re- 
gion centered around r nuc on the disk midplane. Note that the 
time-averaged thickness of this region (~ r nuc ) is larger than 
the instantaneous thickness of the burning region in the case 
of a detonation (Fig. 4). 

Figure 1 1 compares the time-averaged disk properties of 
otherwise identical models with and without nuclear burning 
(COq025_nc and COq000_nc, respectively). The most sig- 
nificant effect of nuclear burning lies in the properties of the 
polar funnel, which is more unbound than in the model with- 
out nuclear burning. This difference also manifests in the en- 
tropy, which is higher inside the funnel for disks with burn- 
ing. A higher entropy indicates excess heating in low temper- 
ature regions, as could result from additional heating above 
the midplane due to the enhanced convection from nuclear 
burning (Fig. 8). An entropy-increasing heating profile is pre- 
cisely the kind necessary to power an energetic outflow, as is 
well known in other astrophysical contexts such as the solar 
wind. 

Nuclear burning does not greatly affect the structure of the 
disk (p > pq\ eq. [36]), although regions with lower density 
are -~50% more extended in radius. The presence of nuclear 



burning also does not alter the fact that surfaces of constant 
entropy and angular momentum are parallel, as is found in 
hydrodynamic simulations of RIAFs (e.g., Stone et al. 1999). 
This indicates that the disk is still marginally stable to con- 
vection, even when the energy input by nuclear heating is in- 
cluded. 

Table 3 gives the time-averaged net accretion rate (M- m ) at 
the inner boundary, total outflow rate, as well as unbound out- 
flow at r = 4/?o in our quiescently-accreting models. Note that 
(A^unb) increases with In non-exploding models with para- 
metric burning rate, the mass loss in unbound material ap- 
proaches a substantial fraction of (M- m ) as ^ — > ^-it. The 
fraction of ash in the unbound outflow also correlates with "J/. 

The decrease of the mass accretion rate with allows the 
disk to survive for a time longer than the characteristic viscous 
time. Model CO_f 1 is an extreme in this sense. The mass-loss 
in unbound flows exceeds the accretion rate by two orders of 
magnitude. Most of the WD material will be ejected over 
^100 orbits, or about 1.5 hours as an unbound flow with a 
velocity ~ 1000 km s" 1 . Converting to physical units, the 
mass loss rate in unbound material can lie in the range 10~ 5 - 
10~ 3 Mq s -1 for the full range of models explored. 



5. DISCUSSION 

5.1. Implications for WD-NS/BH mergers 

An obvious implication of our results is that the C-O torii 
created by WD-NS can undergo a large-scale detonation and 
explosion. This explosive evolution differs drastically from 
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the steady-state model for NuDAFs envisioned by Ml 2, 10 al- 
though the Ml 2 model does provides a qualitative descrip- 
tion of the quiescent mode of accretion found to occur when 
$ < Sl/crit- A similar conclusion would likely apply to the case 
of a hybrid He-C-O WD because the reaction 4 He( 16 O,7) 20 Ne 
has a similar threshold ^ cr i t for detonation (Table 2). Given 
the ubiquity of turbulence generated by the MRI in physical 
accretion disks and relative ease with which the requisite con- 
ditions for a detonation are met (Appendix C), it is possible 
that detonation is a robust outcome of any disk with \1/ > 1 
under the influence of a sufficiently temperature-sensitive nu- 
clear reaction. Our picture of turbulence-generated detona- 
tions in many ways resembles that invoked in 'pulsation de- 
layed detonation' models for Type la SNe (e.g., Khokhlov 
etal. 1997;Niemeyer& Woosley 1997). 11 

Although detonations appear to be a common feature of 
disk evolution when "J > 1, it is challenging to make a defini- 
tive statement about the ubiquity of a detonation because there 
are hints that the outcome could depend on more parameters 
than alone. For instance, our models CO_fl and CO_f2 
with full reaction rates included have relatively minor differ- 
ences in their parameters, yet one explodes while the other 
does not. The precise composition of the WD may also mat- 
ter since $ is proportional to the fuel mass fraction Xf. 

Another important caveat is that our current EOS includes 
only ideal gas pressure, despite the fact that radiation pressure 
is important at small radii and early times in the disk evolu- 
tion. Since radiation pressure acts to reduce the temperature 
at a given pressure, this implies that (1) a deeper potential 
well is required for burning, thus reducing the value of "J; and 
(2) the likelihood that runaway burning will occur is also re- 
duced since the magnitude of temperature fluctuations across 
the burning front are suppressed (Appendix C). It is thus pos- 
sible that radiation pressure could suppress a detonation at 
early times in the torus evolution. However, this should not be 
the case once steady-state accretion is achieved on a timescale 
~ fvisc since then the density of the inner accretion flow is suf- 
ficient for gas pressure to dominate. The radially-decreasing 
density profile achieved on a timescale t > f v j sc is also the most 
conducive to a sustaining an outward-propagating detonation 
(§4.2). This gives us confidence that the threshold for detona- 
tion found in our simulations is not an artifact of our EOS. 

One consequence of a relatively prompt detonation is that a 
large fraction of the WD is unbound, which substantially lim- 
its the mass ultimately accreted by the central NS. This makes 
it less likely that a black hole will be created following a WD- 
NS merger (Paschalidis et al. 201 1). A prompt explosion also 
limits the energy of a relativistic outflow (i.e. a 'jet') orig- 
inating from the vicinity of the compact object (King et al. 
2007) and its resulting high energy or radio emission. Even 
when a strong detonation occurs, however, the existence of a 
remnant disk of bound material (Fig. 6) indicates that a small 
mass ~ 10~ 2 M Q is still available to accrete. 

10 MI2 did show that his solutions were prone to thermal instability due 
to the sensitive temperature dependence of nuclear reactions. Although one 
interpretation of our results is that such instabilities manifest as a detonation, 
the Ml 2 model only accounts for the properties of the mean flow and not for 
the role of turbulence, which plays such an essential role in the detonation 
process in our simulations. 

1 1 How the deflagration ignited near the core of a Chandrasekhar-mass C- 
O WD transitions into a detonation remains a major theoretical uncertainty 
in standard Type la supernova models. Assuming that the initial deflagra- 
tion fails to unbind the WD, delayed detonation models invoke a secondary 
detonation which results from the mixing of partially-burned fuel and ash 
following the subsequent contraction of the WD. 



Our highly idealized initial conditions are necessarily dif- 
ferent from torii that arise in self-consistent mergers. On the 
one hand, the disruption process is not instantaneous (e.g., 
Fryer et al. 1999), thus the angular momentum distribution 
will not be constant. Also, the thermalization of the rotational 
kinetic energy by shocks likely leads to a non-uniform en- 
tropy distribution in the disk. If detonations early in the accre- 
tion phase turn out to be a prevalent outcome, then the initial 
form of the torus will be an important factor to consider when 
making observational predictions. This would also imply that 
the angular momentum transport process must be treated cor- 
rectly. However, we do not expect some of the more general 
results of this paper, such as the existence of a critical value of 
& ~ 1 for detonation or the dependence of quiescent outflows 
on \I>, to depend fundamentally on the details of the initial 
condition. 

In the case of a WD-BH merger, the larger mass of the cen- 
tral compact object results in a lower value of ^ oc M~ l , mak- 
ing it more likely that accretion will instead occur in the qui- 
escent regime (^> < ^crit). If non-explosive burning allows a 
larger fraction of the WD to accrete, then the power of an as- 
sociated high energy counterpart (Fryer et al. 1999) may be 
larger than in the WD-NS case. One limitation on the rate 
of such events is that a binary with a high mass ratio is re- 
quired for unstable mass transfer to occur in the first place 
(q > 0.2-0.5; Paschalidis et al. 2009). The low mass black 
holes (M c < 5Mq) thus required appear to be relatively rare 
among the population of Galactic binaries (e.g., Bailyn et al. 
1998; Fryer et al. 2012). 

In contrast to C-O torii, He torii appear unlikely to explode, 
even under conditions we find conducive to a detonation (Ta- 
bles 1 and 3). Stability results mostly from the weaker tem- 
perature dependence of the triple-a nuclear reaction rate at 
high temperature, which makes the conditions for detonation 
via turbulent mixing less likely to be satisfied (Appendix C). 
The torus density in our simulation is sufficiently low that a 
large fraction of 4 He may reach sufficiently high temperature 
to photodissociate before burning into heavier elements. Our 
conclusion that He torii are unlikely to detonate is at odd with 
the suggestion of Schwab et al. (2012) that such detonations 
might occur in disks created by WD-WD mergers. One possi- 
ble difference in their case could be the importance of degen- 
eracy pressure or the higher densities achieved by compres- 
sion due to the presence of a WD surface. 

Another limitation of applying our current simulations di- 
rectly to physical mergers is that we include only a single re- 
action. In some cases this could affect our conclusions about 
the likelihood of a detonation. Although pure helium torii ap- 
pear to be stable, the products of He burning ( 12 C or 16 0) are 
themselves prone to explosive burning. In C-O torii, on the 
other hand, subsequent reactions (e.g., oxygen burning) will 
generally occur much deeper in the potential well (lower r nuc 
and hence lower 'J; Table 2), such that the first reaction has 
the largest dynamical impact. Nevertheless, including a full 
reaction network could result in a much richer evolutionary 
history, with, for instance, several potential detonations oc- 
curing in succession as matter slowly sinks deeper into the 
potential well and burns to increasing heavier elements. This 
behavior may qualitiatively resemble the late stages of un- 
stable shell burning in massive stars (e.g., Arnett & Meakin 
2011; Quataert & Shiode 2012) or pulsational pair instabili- 
tity supernovae (e.g., Woosley et al. 2007). 
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5.2. Supernova-Like Optical Transients 

M12 proposed that the ejecta from a WD-NS/BH merger 
could produce an optical transient, powered by the radioactive 
decay of 56 Ni synthesized in the accretion disk. He was mo- 
tivated by the recent discovery of several Type I supernovae 
which are dimmer and/or more rapidly evolving than normal 
SNe la or Ib/c (e.g., Li et al. 2003; Jha et al. 2006; Foley 
et al. 2009; Perets et al. 2010; Waldman et al. 2011; Kasli- 
wal et al. 2011). The rapid evolution and low luminosities 
of these events require both lower 56 Ni yield and lower total 
ejecta mass than values characteristic of normal SNe. Some of 
these events occur in locations far outside of their host galax- 
ies (Perets et al. 2010; Kasliwal et al. 2011), possibly consis- 
tent with the location of a WD-NS merger if the NS was given 
a natal 'kick'. 

Based on a steady-state model of accretion following a 
WD/NS-BH merger, M12 predicted an ejecta composed pri- 
marily of unburned C-O or He (depending on the initial WD) 
along with a range of intermediate mass elements and a small 
quantity of 56 Ni. Since our calculations in this paper include 
only a single reaction, we cannot directly predict the Ni yield 
or precisely determine the velocity structure of the ejecta. 
Nevertheless, we can address whether qualitative features of 
the outflows are consistent with those required to produce a 
supernova-like transient. 

In the case of quiescent disk evolution (§4.3), the picture 
that we find is qualitatively similar to that of M12. Nuclear re- 
actions power a quasi-steady outflow with an enhanced mass- 
loss rate over the case without nuclear burning (Fig. 11). If 
* < ^cnt is satisfied for the first reaction activated at large 
radii in the disk, then this condition will also be satisfied 
for subsequent reactions which release at most a compara- 
ble amount of energy but occur deeper within the potential 
well. The total mass of 56 Ni ejected is indeed likely to be 
small Mwi ~ 10~ 3 - 10~ 2 M Q , because mass loss from the outer 
disk cuts off the supply reaching smaller radii where the 56 Ni 
produced. Note, however, that since our simulations show 
that unbound outflows from the disk primarily originate from 
regions interior to the burning radius (Figs. 9, 10), this sug- 
gests that the M12 model may overestimate the fraction of 
unburned fuel in the ejecta. 

On the other hand, if the disk undergoes a global detona- 
tion (§4.2), then this picture is drastically altered. In this case 
the yield of intermediate mass elements and 56 Ni will instead 
depend on how far nuclear reactions proceed behind the det- 
onation front. Since the average density of the outer disk is 
relatively low compared to that of the original WD, the mass 
of 56 Ni synthesized will almost certainly be much less than in 
a normal Type la SNe (cf. Sim et al. 2010). How much Ni is 
produced will depend on the time of detonation fdet relative to 
the viscous time of the torus f v i sc (eq. [8]). If fdet <C fvisc, then 
only a small fraction of the torus mass has spread to small 
radii where the density is sufficiently high for 56 Ni to be pro- 
duced, while if f(j et ~ f v j sc then an larger fraction of the shocked 
WD material will be processed to 56 Ni. Our current simula- 
tions show that the former case (fdet <C fvisc) is more likely, 
but it is still possible that the latter case is more physical if 
radiation pressure indeed suppresses early detonations. 

5.2.1. Connection to SN 2002cx-like Events ? 

One type of SNe with characteristics seemingly compatible 
with those resulting from a WD-NS/BH merger are the events 



prototyped by SN 2002cx and SN 2005bj ("SN 2002cx-like" 
events; Li et al. 2003; Jha et al. 2006; Phillips et al. 2007; 
Foley et al. 2009). This rare class of la SNe are in part distin- 
guished by the following features: (1) high ionization spec- 
trum at maximum light, and the presence of elements lighter 
than the Fe group at all epochs, both indicating that the ejecta 
is extensively mixed; (2) a low peak luminosity compared to 
that expected given its rate of decline (i.e. not following the 
standard 'Phillips relation'), indicating a wide range of syn- 
thesized 56 Ni mass, ranging fromM Ni ~3x 10" 3 -0.2M Q ; (3) 
low expansion velocities ~ 3000-5000 km s" 1 (hence low ki- 
netic energy) compared to a normal SN la (~ 10, 000 km s -1 ); 
(4) permitted Fe II lines and continuum photospheric emis- 
sion at late times, indicating very low velocities < 1 , 000 km 
s -1 for the innermost ejecta; and (5) host galaxies with both 
active star formation (Foley et al. 2009) as well as older stellar 
populations (Foley et al. 2010). (6) An estimated occurence 
rate < 10% of that of normal SNe la (e.g. Phillips et al. 2007). 

Many of the above characteristics are qualitatively consis- 
tent with the expected SN produced by a disk detonation fol- 
lowing a WD-NS or WD-BH merger. Given that a wide range 
in disk densities is expected depending on the mass ratio of 
the binary, one likewise expects a wide range in the mass 
of 56 Ni and intermediate mass elements. Extensive mixing 
also seems plausible given the asymmetric nature of the ex- 
plosion. Low velocity of the inner ejecta is also predicted, 
since the inner ejecta is substantially slowed by the gravita- 
tional field of the central compact object (not present in a 
usual Type la event). Although the expected population of 
host galaxies is difficult to predict with confidence, a mixture 
of both early and late-type galaxies is likely (e.g. Belczyn- 
ski et al. 2002), given the distribution of gravitational wave 
inspiral times. The lack of spectroscopic evidence for un- 
burned 12 C in SN 2002cx-like events is also compatible with 
a disk detonation model, though this observation is in tension 
with alternative models invoking the pure deflagration of a 
Chandrasekhar-mass WD (e.g., Blinnikov et al. 2006). The 
rate of WD-NS mergers in the Milky Way is also estimated 
to be ~ 10~ 4 yr" 1 (based on known population of tight WD- 
NS binaries; O'Shaughnessy & Kim 2010), or approximately 
~ 3% of the SN la rate. One potential problem with the disk 
detonation model is the continuum polarization, which one 
would naively expect to be high given the large asymmetry, 
yet at least in one case (2005hk) was observed to be small 
(~ 0.4%; Chornock et al. 2006). 

We will explore the properties of supernovae from WD- 
NS/BH mergers in greater detail in future work using 
the properties of outflows from the torus including a full 
a-reaction network. These will allow us to better test 
whether WD-NS/BH mergers can be associated with 2002cx- 
like events or other classes of subluminous Type I SNe. 

5.3. Implications for Collapsar Accretion Disks 

Another application of our results is to accretion disks 
formed by the collapse of rotating stars, as in collapsar models 
for gamma-ray bursts (GRBs) (MacFadyen & Woosley 1999; 
Lindner et al. 2012). Depending on the angular momentum 
profile of the star (Woosley & Heger 2006), a significant frac- 
tion of the collapsing envelope may circularize at sufficiently 
large radii > few x 10 8 cm, where nuclear reactions such as 
4 He( 16 O,7) 20 Ne are important. Given the relatively low mass 
of the black hole formed by the collapse of a Wolf-Rayet star, 
one expects ^ <~ ^ cr i t in collapsar disks, such that nuclear 
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burning could indeed impact the dynamics of the accretion 
flow. If the effects of nuclear burning produce large ampli- 
tude fluctuations in the central accretion rate, then the viscous 
or dynamical timescale near the burning radius (~ seconds 
for r nuc ~ few x 10 8 cm) could be imprinted on the variabil- 
ity of GRB emission. This could explain, for instance, why 
the power density spectra of GRB light curves peaks at a fre- 
quency ~ Hz (Beloborodov et al. 2000) which is much lower 
than the ~ kHz variability characteristic of the innermost sta- 
ble orbit. 

Quiescent outflows, or periodic episodes of runaway burn- 
ing, could also provide a source of 56 Ni, as is required to 
power the light curves of GRB supernovae. Bodenheimer & 
Woosley (1983) first showed that a sufficiently rapidly rotat- 
ing star could undergo explosive oxygen burning upon col- 
lapse, due to heating caused by the centrifugal 'hang up' of in- 
falling material (cf. MacFadyen & Woosley 1999). How much 
56 Ni is synthesized from such explosive burning of stellar fuel, 
versus that produced within much hotter outflows from the in- 
ner torus (MacFadyen & Woosley 1999; Metzger et al. 2008; 
Milosavljevic et al. 2012; Surman et al. 2011), has yet to be 
quantified. 

Nuclear burning may equally be important in accretion fol- 
lowing the merger of a He star with a NS or BH (Fryer & 
Woosley 1998). Thone et al. (2011) proposed that a He star- 
NS merger was responsible for the highly unusual "Christ- 
mas" gamma-ray burst 101225 A, which exhibited exception- 
ally long gamma-ray and thermal X-ray emission. In order to 
explain the peak luminosity of the optical 'bump' following 
this event as supernova-like emission, the ejected mass of 56 Ni 
must have been small < 0. 1M Q (if one adopts the distance ad- 
vocated by Thone et al. 201 1). Future work on NuDAFs will 
better address the accretion efficiency and 56 Ni ejected due to 
rapid He accretion, thereby allowing us to assess whether the 
high energy and thermal optical emission from GRB 101225 A 
was indeed compatible with a He star-NS merger. 

6. SUMMARY 

This paper explores the effect of nuclear reactions on the 
evolution of RIAFs, in the context of merging WDs and 
NSs or BHs. Two-dimensional hydrodynamic simulations 
with FLASH3.2 are used to systematically characterize the 
properties of these disks. Our main findings are the following: 

1 . - The effect of nuclear burning on RIAFs is controlled by 
the ratio "J of the nuclear energy to the enthalpy at the radius 
r nuc where most of the fuel is consumed (eqns. [37] and [39]). 
The qualitative behavior of the system depends on the value 
of relative to a critical value ^> c ^ t ~ 1 which separates 
quiescent burning from large-scale detonations. The exact 
value of this critical parameter is sensitive to details about 
the disk, such as the thermal content or the amplitude of the 
viscous stress (Table 3). 

2. - For disks that include cooling and/or have \I> > ^crit, 
a detonation can be triggered by hot spots formed near the 
burning front (Fig. 3). These temperature enhancements 
are generated by turbulent mixing of cold fuel and hot 
ash (Fig. 4), which produce an induction time-gradient 
consistent with that required to generate a detonation via the 
Zel'dovich mechanism (Fig. 4; Appendix C). For low values 
of detonations remain localized, and provide at most a 
slight enhancement to the mass ejection in quiescent disks. 



Increasing gives rise to detonations of increasing power. 

3. - Exploding disks ("f ^> & C nt) can easily achieve expansion 
velocities in excess of 1,000 km s" 1 (Table 3). A remnant 
disk, with a mass typically a few percent of the initial WD 
mass, is left behind (Fig. 6). 

4. - Non-exploding disks with nuclear reactions can generate 
unbound outflows along a funnel next to the rotation axis 
(Fig. 9). This material is composed primarily of nuclear ash. 
The expansion velocity of these outflows can also achieve 
~ 1000 km s -1 (Fig. 10), with mass outflow rates spanning 
the range 10" 5 - 10" 3 M s" 1 (Table 3). 

5. - The energy deposition by nuclear reactions locally 
enhances the turbulent kinetic energy, decreasing the mass 
accretion rate with increasing (Fig. 8). Given that a 
significant fraction of the outflowing material is bound and 
hence returns to the disk (Fig. 9), this choking of accretion 
can prolong the lifetime of the disk well beyond a few 
viscous times. A significant amount of burnt material can 
then accumulate at large distances from the disk. 

6. - Outflows from the disk, generated either via quiescient 
disk winds or a large-scale detonation, may give rise to 
an optical supernova powered by the radioactive decay of 
56 Ni. In the case of a disk detonation, some of the properties 
of the predicted transient are consistent with those of the 
observed class of unusual Type la SNe defined by SN 2002cx. 

This study has focused primarily on characterizing the 
range of outcomes obtained in different regions of parameter 
space, and on identifying the main parameter dependencies. 
Given the number of approximations made, our results cannot 
be directly translated into observational predictions, nor do 
they constitute definitive statements about the likely outcome 
of a realistic WD/NS or WD/BH merger. Aside from includ- 
ing a realistic equation of state and a full nuclear reaction net- 
work, a convergence study on the circum-torus medium needs 
to be performed in order to reliably predict expansion veloc- 
ities and thermodynamic properties of the ejecta. Also, fully 
three-dimensional simulations are eventually desirable, since 
the formation and ignition of hot spots is likely to occur at a 
position which is well-localized in azimuth, an effect clearly 
not captured by an axisymmetric calculation. 

Even if NuDAFs in nature turn out not to explode, the fact 
that realistic burning rates yield ^ > 1 (Table 2) implies that 
the structure and evolution of these disks will indeed be rich 
due to the multiple elements to be burned and the dynamical 
importance of each of these reactions. For example, despite 
the fact that model CO_fl [full 12 C( 12 C,7) 24 Mg rate] does not 
achieve explosion, the mass loss rate in the unbound material 
is ~ 100 times the net accretion rate at the inner boundary. A 
significant fraction of the disk will be ejected as an unbound 
flow at high velocities over a timescale ~ 100f or b ~ few hr. 
The observational signature of such an event may constitute a 
unique type of transient, which is distinct from that produced 
in the case of disk detonation or as yet observed. 

A companion paper will continue the study of these disks 
using a realistic equation of state and a full nuclear reaction 
network. 
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APPENDIX 

A. UNIT CONVERSION AND BURNING RATES 

For completeness, this Appendix show explicitly the conversion from physical to code units and the functional forms of the 
reactions used. 

The code takes the circularization radius Rq, the orbital velocity at this radius (GM c /Rq) 1 I 2 , and the maximum initial torus 
density 

f< ln5 /Pmax\ / M WD \ /I0"cm\ 3 _ 3 

as the basic set of units. In equation (Al), (p max / p) is the ratio of the maximum to average density in the torus. This quantity is 
obtained by computing the total torus mass using the density profile in equation (26), being a function of the adiabatic index and 
distortion parameter only. For 7 = 5/3 and d = {1.2, 1.5,3}, the result is yO max /p = {0.536,0.176,2.47 x 10~ 2 }. The temperature 
corresponding to p/p = GM c /Rq is 



r -= 2x n^Ju^A^j K ' (A2) 

where fj, is given by equation (7). The mean molecular weight for temperature conversion purposes is calculated using the initial 
composition of the disk, once per model. Given the pressure and density in code units, we use equations (Al) and (A2) to obtain 
physical density and temperature 

The specific energy generation rates in physical units are converted to code units through division by 

1 fGM c \ V2 . . M c ^ 3/2 



M c \ /10 9 ' 3 



cm 



1.4 x 10 
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The rate of change of mass fractions is converted to code units through multiplication by 

Rl \ 1/2 / R \ 3/2 /1.4M Q \ 1/2 

Analytic nuclear burning rates are taken from Caughlan & Fowler (1988) 12 . The bulk of our study makes use of the 
12 C( 12 C,7) 24 Mg reaction, which is the most energetic among the a reactions involving carbon, oxygen, and helium. The specific 
energy generation rate is 



T 5/6 

enuc,cl2 = 3.96xl0 43 pi^C^2 
^9 



ex P y-^TJT ~ 2 - 12 x W~% 3 J er g g~' s "' > ( A 5) 



where p\ is the density in g cm , Tg = T/(10 K), Xc is the mass fraction of C, and 



T A g= . (A6) 

l+0.0396r 9 



The rate of change of Xq is given by 
where mc = 12m„ is the mass of a carbon nucleus, and gi2 = 13.933 MeV is the energy liberated in the reaction 



^C,cl2 =-77— 2nuc,12, (A7) 



12 http : / /www . phy . ornl . gov/ astrophysics /data /cf 8 8/ 
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To gain an analytic understanding of the effect of nuclear burning, we also use an approximate power-law form of the 
12 C( 12 C,7) 24 Mg reaction. Following a standard procedure (e.g., Kippenhahn & Weigert 1994), we expand equation (A5) in a 
Taylor series in temperature around the point where nuclear burning is equal to the dynamical time, and obtain 

2nuc.pi = 3.06 x 10 6 X 2 p\ r 9 29 erg g" 1 s" 1 . (A8) 

The change in the carbon mass fraction is obtained from equation (A7) by replacing the energy generation rate. Good agreement 
between the two formulations is found in the temperature range [0.6, 1.2] x 10 9 K. At higher temperatures, the full rate has an 
increasingly weaker dependence on temperature relative to the power-law approximation. 
To consider the case of helium white dwarfs, we also include the triple-alpha reaction. The energy generation rate is 

<2nuc,3a = 5.07 x 10*X^p] T 9 ~ 3 exp 

+2.45 x 10 9 r 2 Ap\ r 9 - 3/2 exp , (A9) 

where Xu e is the mass fraction in helium nuclei. The coefficient r28 is set to 1 /10. The rate of change of the helium mass fraction 
is 

^He = - ^— 2nuc,3a, (A10) 

mc 

where Qj, a = 7.275 MeV is the energy released in the reaction. 
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B. VERIFICATION TESTS 

This Appendix describes a series of tests conducted on the implementation of the NuDAF setup in FLASH3.2 (§3.3). 

We first tested the degree to which the torus can remain in steady-state in the absence of viscous source terms. To this end, we 
use a torus with d = 1.125, and an isothermal atmosphere with p. dt = 10~ 6 /9 max at r at = QARq (eq. [28]). The computational domain 
extends from r- m = QARq to r out = 4Rq, and the resolution is N r = 64, with 90 uniformly spaced cells covering the full range of polar 
angles. Numerical diffusion causes a very small amount of mass to peel off from the torus edges and accrete through the inner 
boundary. At our baseline resolution (§3.3), the total angular momentum in the computational domain is conserved to within 
a few parts in 10~ 6 over 10 orbits. The mass and angular momentum inside torii (defined as the material inside an iso-density 
surface at 10~ 3 /9 max ), are conserved to within a few parts in 10~ 4 over 10 orbits. Over the same time interval, the maximum torus 
density decreases by slightly more than 1%. 

The inner boundary condition (eq. [33]) was tested by evolving a torus similar to that used in the previous test, but now with 
neither angular momentum nor viscosity. Under such conditions, the system proceeds to free-fall towards the gravitating mass. 
For an inner boundary with diameter comparable to the torus thickness (n n = 0.4/?o), the torus material is cleanly absorbed by 
the inner radial boundary, with no transients being generated. Using a smaller inner radius (r m = 0.04/?o) causes some material to 
collide with the symmetry axis, but the bulk of the torus is still cleanly accreted, with no discernible feedback from the boundary. 

To test the accuracy of our viscous diffusion operator, we multiply equation (22) by £ z , make use of equations (20) and (29), 
and integrate over volume, to obtain 

J^L 2 + /d 3 *F r V4 = 0, (Bl) 

where 

L 2 = J d 3 x (iplfj , (B2) 

and the integral is carried out over the full computational domain. For equation (Bl) to be satisfied, the advection of l z and the 
flux formulation of the viscous operator must be treated correctly (T. Heinemann, private communication). Figure 12 shows the 
relative difference between the two terms in equation (Bl), averaged over the first orbit, for a test model evolved with different 
spatial resolutions. The torus has d = 1.125 and a viscosity proportional to density (eq. [15]), with vq = 0.01. Values lie between 
10~ 3 and 10~ 2 depending on resolution. A number of operations are involved in computing this equation, which degrade its 
accuracy, thus we consider the results as indicative of our correct treatment of viscous diffusion within the limitations of the 
hydrodynamic method. 

As an additional test of angular momentum transport, we compared the result of using the flux-conservative formulation for 
angular momentum transport (eqns. [29] and [30]) with the result of directly computing the divergence of the viscous stress tensor 

( V • Y)+ = \ |- {r 3 T I<p ) + -1- ^ (sin 2 6 T g<p ) (B3) 
ror v ' sin 8 oO v ' 

as a 3-point finite difference operator. The result differs at the edges and surroundings of the torus, with the flux formulation 
maintaining a sharper torus surface (in £ : ) for a longer time. Inside the torus, the evolution is nearly identical within stochastic 
fluctuations induced by convection. 
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FIG. 12. — Relative difference between the two terms appearing in equation (Bl), normalized relative to the second term, and averaged between 20% and 80% 
of the first orbit. The torus has d = 1.125, inner and outer radii r m = OARq and r ou t = 4Rq, respectively, and a viscosity given by equation (15) with i>q = 0.01. 
The quantities N r and Ng denote the number of grid cells per decade in radius and in the polar direction, respectively. The error bars denote the root-mean-square 
fluctuation of the difference, which is due to fluctuations in the time derivative of equation (B2), computed as a time-centered finite difference with time step 
equal to 1 /1000 of an orbital period. 



C. CONDITIONS FOR DETONATION VIA TURBULENT MIXING 

In this Appendix we evaluate whether the turbulent mixing of ash and fuel can generate the conditions for a detonation 
in NuDAFs. The ZeFdovich criterion for spontaneous initiation of a detonation is that disparate portions of the burning re- 
gion be separated by a distance such that the difference in the timescale for nuclear heating implies a supersonic phase ve- 
locity (Zel'dovich et al. 1970; Blinnikov & Khokhlov 1987; Woosley 1990). Quantitatively, this condition is expressed as 
(e.g., Niemeyer & Woosley 1997; Bell et al. 2004) 

iV^df 1 > c„ (CI) 
where c s is the adiabatic sound speed and fj n d is the induction time, which is loosely defined as the time required to 'run away' 
to high temperatures via burning at constant pressure. Since fi„d oc c 2 / Q auc is a rapidly decreasing function of temperature, a 
shallow temperature gradient is required across the burning region for triggering a detonation. Such a condition may be satisfied 
in NuDAFs due to the effects of turbulent mixing. 

We focus on the burning of fuel with an initial mass fraction XfQ. The steady-state structure of the accretion flow is characterized 
by two radially-separated regions: (1) unburned fuel with temperature 7f and density pf, and (2) burned ash with temperature T. d 
and density p a . The burning front separating the upstream (region 1) from the downstream (region 2) is centered about the radius 
r = r nac (eq. [39]) and has a width Ar nuc . The radius r nuc is approximately determined by equality of the fuel consumption time 
W = Xf/Xf and the dynamical timescale fd yn = (r^/GM) 1 ' 2 . 

If burning occurs at constant pressure, then conservation of mass and enthalpy determines the change in density and temperature 
across the burning front (e.g., Khokhlov et al. 1997): 



1 



Pf T, 1 + * 

Pa 1 



(/gas > Aad) , 



Pf 1+\I/ 



; r f = r, (P^^Pgas), 



where 



<f 



(C2) 
(C3) 

(C4) 



is the ratio of the specific nuclear energy released per reaction (eq. [19]) to the enthalpy of the fuel c 2 f /(j- 1) (see also equa- 
tion [37]). We have separated cases corresponding to whether gas or radiation pressure dominates. 

Due to the turbulent nature of accretion and RT instabilities, the burning front at r = r nuc is not completely smooth (e.g., Fig. 3). 
Instead, parcels of burned ash occasionally mix with the colder fuel upstream. If mixing produces conditions such that equation 
(CI) is satisfied over the length of the eddy L e dd, then a detonation may be triggered. If a fraction / a of ash is mixed with a 
fraction 1 — /„ of fuel (resulting in a reactant mass fraction X m = [1 -/ a ]Xfo), then the density p m and temperature T m of the 
resulting mixture are also determined by conservation of mass and enthalpy: 



Pm _ Tf_ _ 1 

Pm 1 



(Pens 3> Prad) , 



Pf l+/a# 



i T m = Tf (P ra d> J Pgas), 



(C5) 

(C6) 
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The length of an eddy L e dd which is capable of efficient mixing must obey Ar nuc < L e dd < H = r nuc /2 since it must both fit 
inside the disk midplane and it must sample both fuel and ash over the width of the mixing region Ar nuc . The condition for 
detonation in equation (CI) can thus approximately be written as 



|Vf ind | 



^edd 



or 



^ind(^m; Pirn -^rn) 

Pm,T m ) 



>C S 



> 



^orb(y nuc) 



(C7) 



(C8) 



where we have assumed that the sound speed of the mixture is similar to that in the disk midplane, and we have used the condition 
of vertical hydrostatic equilibrium c s =H/t OI b. 

Assuming a burning rate of the general form Q nuc oc e nuc X q p q ~ x T^ (eq. [A8]), one has that f; n d oc p~te~ 1 )X~ 9 T~P /[$> /3] for 
burning at constant pressure, where the factor of /? in the denominator accounts for the fact that runaway to high temperatures 
occurs on a timescale which is shorter than the heating timescale estimated using the initial heating rate when gas pressure 



dominates (/3 ~ j3 when P gas > P nd ; /3 ~ 1 when P lad > P gas ). Since f nuc = X/X oc p-^-^X'^'^T' 13 , we also have f, 
Using the relation t lwc (Xf, pf, Tf) = td yn that determines r nuc , one can rewrite equation (C8) as 



ind • 



Ledd 1 
"7T > 27T/3* 



Pm,T m ) 
tind(Xf,Pf,Tf) 



(C9) 



or 



^edd 

H 




X m 



(l+/a*) 9 , 



(Pgas/Prad) > 1 
(W-Pgas) > 1 ' 



(CIO) 



where we have used the results from equation (C6). 

Equation (CIO) shows that in a gas pressure-dominated flow, even small eddies with L e dd ~ 0.1//, such as those producing 
detonations in our simulations (Fig. 3), require only a moderate amount of entrained fuel f a <C 1 to satisfy the condition for a 
detonation given a temperature-sensitive reaction (e.g., f3 = 29, q = 2; eq. [A8]). The mixed fraction of upstream ash is indeed 
small (X a ~ few %) in the regions which form hot spots and detonations in our simulations (see Fig. 4). The sensitive dependence 
of the bracketed quantity on the RHS on "J/ (almost exponential) may in part also explain why the formation of a detonation in 
our simulations depends sensitively on this value. On the other hand, in a radiation pressure dominated disk the condition (CIO) 
is more challenging to satisfy, as the lack of temperature fluctuations cause mixing to be detrimental given the decrease in density 
(eq. [C6]). 

Producing a detonation is harder than satisfying condition (CIO), since mixing must occur before complete burning. Another 
important parameter is thus the ratio of the characteristic eddy turnover timescale f e dd ~ /^ddAedd = -M -1 L el \&/c s , which sets the 
timescale for mixing, to the burning timescale of the mixture ? nuc : 



fedd \£edd/?ind 



(Cll) 



where A4 = Vedd/c s < 1 is the turbulent Mach number. Provided that the condition for detonation (eq. [C8]) is marginally 
achieved, equation (Cll) shows that efficient mixing requires a combination of vigorous turbulence and/or strong energy release. 
Dominance of radiation pressure makes mixing more inefficient, which is helpful in satisfying equation (CIO), as a higher mixing 
fraction results in a larger required eddy size. 



REFERENCES 



Abramowicz, M. A., Lasota, J.-R, & Igumenshchev, I. V. 2000, MNRAS, 
314, 775 

Arnett, W. D., & Meakin, C. 201 1, ApJ, 733, 78 

Bailyn, C. D., Jain, R. K., Coppi, P., & Orosz, J. A. 1998, ApJ, 499, 367 
Balbus, S. A., & Hawley, J. F. 1998, Reviews of Modern Physics, 70, 1 
Belczynski, K., Bulik, T„ & Rudak, B. 2002, ApJ, 571, 394 
Bell, J. B., Day, M. S., Rendleman, C. A., Woosley, S. E„ & Zingale, M. 

2004, ApJ, 608, 883 
Beloborodov, A. M., Stern, B. E„ & Svensson, R. 2000, ApJ, 535, 158 
Blandfbrd, R. D., & Begelman, M. C. 1999, MNRAS, 303, LI 
Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 
Blinnikov, S. I., & Khokhlov, A. M. 1987, Soviet Astronomy Letters, 13, 

364 

Blinnikov, S. I., Rbpke, F. K., Sorokina, E. I., Gieseler, M., Reinecke, M., 
Travaglio, C, Hillebrandt, W., & Stritzinger, M. 2006, A&A, 453, 229 
Bodenheimer, P., & Woosley, S. E. 1983, ApJ, 269, 281 



Caughlan, G. R., & Fowler, W. A. 1988, Atomic Data and Nuclear Data 

Tables, 40, 283 
Chen, W.-X., & Beloborodov, A. M. 2007, ApJ, 657, 383 
Chevalier, R. A. 1993, ApJ, 411, L33 

Chornock, R., Filippenko, A. V., Branch, D., Foley, R. J., Jha, S., & Li, W. 

2006, PASP, 118,722 
Colella, P., & Woodward, P. R. 1984, JCP, 54, 174 
Davis, S. W., Stone, J. M., & Pessah, M. E. 2010, ApJ, 713, 52 
Di Matteo, T., Perna, R., & Narayan, R. 2002, ApJ, 579, 706 
Dubey, A., Antypas, K., Ganapathy, M. K., Reid, L. B., Riley, K., Sheeler, 

D., Siegel, A., & Weide, K. 2009, J. Par. Comp., 35, 512 
Eggleton, P. P. 1983, ApJ, 268, 368 
Fernandez, R. 2012, ApJ, 749, 142 

Foley, R. J., Chornock, R., Filippenko, A. V., Ganeshalingam, M., Kirshner, 
R. P., Li, W., Cenko, S. B., Challis, P. J., Friedman, A. S., Modjaz, M., 
Silverman, J. M., & Wood-Vasey, W. M. 2009, AJ, 138, 376 



NUCLEAR DOMINATED ACCRETION FLOWS IN 2D 



23 



Foley, R. J., Rest, A., Stritzinger, M., Pignata, G., Anderson, J. P., Hamuy, 
M, Morrell, N. I., Phillips, M. M., & Salgado, F. 2010, AJ, 140, 1321 

Fryer, C. L., Belczynski, K., Wiktorowicz, G., Dominik, M., Kalogera, V., & 
Holz, D. E. 2012, ApJ, 749, 91 

Fryer, C. L., & Woosley, S. E. 1998, ApJ, 502, L9 

Fryer, C. L., Woosley, S. E., Herant, M., & Davies, M. B. 1999, ApJ, 520, 
650 

Hawley, J. F, & Balbus, S. A. 2002, ApJ, 573, 738 
Hawley, J. F., Balbus, S. A., & Stone, J. M. 2001, ApJ, 554, L49 
Hawley, J. F, Gammie, C. F, & Balbus, S. A. 1995, ApJ, 440, 742 
Hirose, S., Krolik, J. H., & Stone, J. M. 2006, ApJ, 640, 901 
Igumenshchev, I. V., & Abramowicz, M. A. 1999, MNRAS, 303, 309 
Igumenshchev, I. V., & Abramowicz, M. A. 2000, ApJS, 130, 463 
Igumenshchev, I. V., Abramowicz, M. A., & Narayan, R. 2000, ApJ, 537, 
L27 

Igumenshchev, I. V., Narayan, R., & Abramowicz, M. A. 2003, ApJ, 592, 
1042 

Itoh, N, Hayashi, H., Nishikawa, A., & Kohyama, Y. 1996, ApJS, 102, 411 
Jha, S., Branch, D., Chornock, R., Foley, R. J., Li, W., Swift, B. J., Casebeer, 

D., & Filippenko, A. V. 2006, AJ, 132, 189 
Kasliwal, M. M., etal. 2011, arXiv:1111.6109 
Khokhlov, A. M., Oran, E. S., & Wheeler, J. C. 1997, ApJ, 478, 678 
King, A., Olsson, E., & Davies, M. B. 2007, MNRAS, 374, L34 
Kippenhahn, R., & Weigert, A. 1994, Stellar Structure and Evolution, 2nd 

edn. (New York: Springer) 
Landau, L. D., & Lifshitz, E. M. 1987, Fluid Mechanics, 2nd edn. (Oxford: 

Butterworfh-Heinemann) 
Larson, R. B. 1984, MNRAS, 206, 197 
Li, W., et al. 2003, PASP, 1 15, 453 

Lindner, C. C., Milosavljevic, M., Couch, S. M., & Kumar, P. 2010, ApJ, 
713, 800 

Lindner, C. C., Milosavljevic, M., Shen, R., & Kumar, P. 2012, ApJ, 750, 
163 

Lisewski, A. M., Hillebrandt, W., Woosley, S. E., Niemeyer, J. C, & 

Kerstein, A. R. 2000, ApJ, 537, 405 
MacFadyen, A. I., & Woosley, S. E. 1999, ApJ, 524, 262 
Machida, M., Matsumoto, R., & Mineshige, S. 2001, PASJ, 53, LI 
McKinney, J. C, Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 

423, 3083 

Metzger, B. D. 2012, MNRAS, 419, 827 

Metzger, B. D., Thompson, T. A., & Quataert, E. 2008, ApJ, 676, 1 130 
Milosavljevic, M., Lindner, C. C., Shen, R., & Kumar, P. 2012, ApJ, 744, 
103 



Narayan, R., Sadowski, A., Penna, R. R, & Kulkarni, A. K. 2012, MNRAS, 

submitted, arXiv: 1206. 1213 
Narayan, R., & Yi, I. 1994, ApJ, 428, L13 
— . 1995, ApJ, 444, 231 
Nauenberg, M. 1972, ApJ, 175, 417 
Niemeyer, J. C, & Woosley, S. E. 1997, ApJ, 475, 740 
Nomoto, K„ & Kondo, Y. 1991, ApJ, 367, L19 
O'Shaughnessy, R., & Kim, C. 2010, ApJ, 715, 230 
Pang, B., Pen, U.-L., Matzner, C. D., Green, S. R., & Liebendorfer, M. 

2011, MNRAS, 415, 1228 
Papaloizou, J. C. B., & Pringle, J. E. 1984, MNRAS, 208, 721 
Paschalidis, V., Liu, Y. T., Etienne, Z., & Shapiro, S. L. 2011, Phys. Rev. D, 

84, 104032 

Paschalidis, V.. MacLeod, M., Baumgarte, T. W., & Shapiro, S. L. 2009, 

Phys. Rev. D, 80, 024006 
Pen, U.-L., Matzner, C. D., & Wong, S. 2003, ApJ, 596, L207 
Perets, H. B., et al. 2010, Nature, 465, 322 
Phillips, M. M., et al. 2007, PASP, 119, 360 
Plewa, T., & Muller, E. 1999, A& A, 342, 179 
Quataert, E., & Gruzinov, A. 2000, ApJ, 539, 809 
Quataert, E., & Shiode, J. 2012, MNRAS, 423, L92 
Schwab, J., Shen, K. J., Quataert, E., Dan, M., & Rosswog, S. 2012, 

MNRAS, submitted, arXiv:1207.0512 
Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 
Shi, J., Krolik, J. H., & Hirose, S. 2010, ApJ, 708, 1716 
Siess, L., & Livio, M. 1999, MNRAS, 304, 925 

Sim, S. A., Ropke, F. K., Hillebrandt, W., Kromer, M., Pakmor, R., Fink, 

M., Ruiter, A. J., & Seitenzahl, I. R. 2010, ApJ, 714, L52 
Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753 
Stone, J. M., & Pringle, J. E. 2001, MNRAS, 322, 461 
Stone, J. M., Pringle, J. E., & Begelman, M. C. 1999, MNRAS, 310, 1002 
Surman, R., McLaughlin, G. C, & Sabbatino, N. 2011, ApJ, 743, 155 
Thone, C. C, et al. 201 1, Nature, 480, 72 

Waldman, R., Sauer, D., Livne, E., Perets, H., Glasner, A., Mazzali, P., 

Truran, J. W., & Gal- Yam, A. 201 1, ApJ, 738, 21 
Wickramasinghe, D. T., & Ferrario, L. 2000, PASP, 112, 873 
Woosley, S. E. 1990, in Supernovae, ed. A. G. Petschek, 182-212 
Woosley, S. E., Blinnikov, S., & Heger, A. 2007, Nature, 450, 390 
Woosley, S. E., & Heger, A. 2006, ApJ, 637, 914 
Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Reviews of Modern 

Physics, 74, 1015 
Yuan, F., Wu, M., & Bu, D. 2012, ApJ, submitted, arXiv: 1206.4157 
Zel'dovich, Y. B., Librovich, V. B., Makhviladze, G. M., & Sivashinskil, 

G. I. 1970, J. Appl. Mech. & Tech. Phys., 11, 264 



